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ABSTRACT 

This  paper  investigates  some  of  the  difficulties 
commonly  encountered  in  the  application  of  pattern  recog- 
nition techniques  to  the  identification  of  non-cooperative 
signals.   Practical  problems  of  feature  selection  and  analy- 
sis are  explored  for  a  data  base  of  actual  bauded  signal 
transition  times.   Analysis  of  a  set  of  features  based  on 
the  raster  display  of  these  transition  times  suggests  that 
the  use  of  fast  Fourier  transform  techniques  may  yield 
useful  identifying  parameters  not  exploitable  by  the  current 
process  of  visual  comparison  of  signal  raster  patterns. 
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I.   INTRODUCTION 

The  manual  measurement  and  analysis  of  signal  parameters 
leading  to  identification  of  the  signal  source  can  require 
more  time  and  manpower  than  results  justify.   Such  a  process 
may  be  made  less  costly  by  the  application  of  computer  tech- 
nology to  measure  the  traditional  parameters  in  a  manner 
relatively  free  of  individual  bias  and  to  aid  in  the  analy- 
sis and  identification  procedures.   An  additional  result  of 
the  computer  application  is  the  capability  to  calculate 
parameters  unmeasurable  by  manual  techniques.   This  thesis 
investigates  the  process  by  which  such  parameters  might  be 
obtained  and  used  in  an  existing  pattern  recognition  scheme. 
The  discussion  progresses  from  the  general  principles  of 
Bayesian  classification  to  the  problem  of  finding  features 
which  are  useful  in  classifying  signals  from  a  specific 
data  base. 

Section  II  of  this  thesis  deals  with  the  theory  of 
pattern  recognition  which  underlies  most  automated  classi- 
fication schemes.   Section  III  discusses  some  of  the  limi- 
tation imposed  on  this  theory  in  a  practical  signal  identi- 
fication problem.   Section  IV  is  devoted  to  the  problems 
which  arise  in  the  attempt  to  choose  a  set  of  features  useful 
and  meaningful  to  the  identification  process.   Section  V 
addresses  the  particular  problems  of  choosing  a  set  of 
useful  features  based  on  the  raster  scan  display  of  bauded 
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signal  transition  times.   Included  in  this  section  are 
the  motivation  for  the  initial  choice  of  parameters  and 
the  techniques  used  for  both  measurement  and  analysis.   The 
conclusions  reached  by  this  investigator,  as  well  as  those 
areas  deserving  further  study,  are  summarized  in  Section  VI 
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II.   BAYESIAN  PATTERN  CLASSIFICATION 

The  process  of  mathematical  pattern  recognition  forma- 
lizes the  methods  by  which  the  measured  parameters  or  fea- 
tures of  a  group  of  samples  may  be  used  to  identify  those 
samples  as  belonging  to  one  of  several  general  classes. 
Such  a  process  may  be  considered  to  be  composed  of  two 
parts.   The  first  part  is  a  learning  procedure  in  which  the 
general  "pattern"  of  each  of  the  various  classes  is  des- 
cribed in  terms  of  the  similarities  and  differences  in  the 
sample  features.   The  second  part  is  a  decision  process  in 
which  a  feature-based  algorithm  is  developed  from  the  pattern 
description  and  applied  to  the  identification  or  classifi- 
cation of  nev;  samples.   Both  parts  of  the  process  are  based 
on  the  theory  of  statistical  inference  and  lead  to  a  pro- 
babilistic interpretation  of  the  class  membership.   Closely 
related  to  such  an  interpretation  is  the  Bayesian  model  of 
decision  theory. 

In  the  Bayesian  model,  the  first  step  of  a  pattern  recog- 
nition scheme,  the  learning  process,  is  that  technique  by 
which  a-priori  class-conditional  feature  probabilities  are 

determined.   The  K . -individual  d-dimensional  feature  vectors 

1 

x,   /  x.o   r  •  •  •  /  X-.    of  the  m  distinct  classes  C.  are 
used  to  estimate  m  class  conditional  probability  density 
functions  on  the  feature  space 
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P(x|ci)  :  p(x|x{i)  •••  zi^)  ("-I) 

i 


Knowing  these  distributions  for  each  class  and  the  overall 
class  probabilities,  p(C),  one  may  obtain  the  feature 
conditional  class  probabilities: 


p(C.)p(x|c.) 

p(C.  |x)  =  1    ~  1- 

r   l1-        p  (x) 

p(C.)p(x|C.) 

=  i_ i (II-2) 

m 

£  p(C.)p(x|C.) 

i=l    X       X 

Assuming  an  equal  cost  of  misclassification ,  the  Bayesian 
decision  rule  chooses  that  class  v;hich  maximizes  the  feature 
conditional  probabilities  of  Equation  II--2.   That  is,  if 

p(C± |x)  >  p(C. |x)     for  all  j  (II-3) 

the  sample  characterized  by  the  feature  vector  x  is  assigned 
to  class  C . . 

The  above  technique  generates  a  series  of  decision 
regions  in  the  d-dimensional  feature  space.   Once  the  a-priori 
probabilities  are  determined,  the  class  selection  is  simply 
an  m-valued  function  of  the  feature  vector.   This  is  essen- 
tially the  goal  of  the  pattern  recognition  process. 

The  probabilistic  Bayesian  decision  process  discussed 
above  is  mathematically  elegant  and  optimal  in  the  sense 
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that  if  the  a-priori  probabilities  are  known,  it  minimizes 
the  expected  risk  of  misclassif ication .   Unfortunately, 
its  implementation  suffers  major  difficulties  in  both  the 
learning  and  the  decision  phases  of  a  statistical  pattern 
recognition  process. 

A.   PATTERN  LEARNING 

In  the  learning  portion  of  a  pattern  recognition  problem, 
one  normally  deals  with  a  finite  number  of  samples  which 
constitute  a  training  set.   These  samples  are  the  basis  from 
which  the  investigator  is  obliged  to  construct  the  necessary 
a-priori  class-conditional  probability  distribution  functions. 
Thus  the  well-defined  probability  structure  assumed  by 
Bayes '  decision  rule  is  in  reality  a  statistical  inference 
based  on  the  data  obtained  from  the  training  set  and  the 
investigator's  best  guess  of  the  functional  form  of  this 
data.   Dependent  on  how  much  of  the  general  form  of  the 
probability  density  function  is  presupposed,  the  estimation 
techniques  are  characterized  as  parametric  or  non-parametric. 

Typically,  parametric  estimation  assumes  that  the  features 
of  the  training  samples  follow  one  of  the  well-known  proba- 
bility distributions.   The  samples  of  the  training  set  are 
used  to  estimate  only  a  few  essential  parameters  of  this 
distribution.   For  example,  one  might  use  the  training  sam- 
ples to  estimate  the  means  and  variances  for  an  assumed 
multivariate  normal  distribution.   The  non-parametric  tech- 
nique on  the  other  hand  seeks  to  estimate  directly  the 
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probability  density  function  as  a  superposition  of  "potential 
function"  terms  arising  from  each  of  the  training  set  samples. 
The  accuracy  of  the  inferred  density  function  is  strongly 
dependent  on  the  type  of  estimation  technique  used,  and  may 
be  severely  limited  by  the  size  of  the  training  set.   The 
applicability  and  limitations  of  both  non-parametric  and 
parametric  estimation  are  discussed  further  in  Appendix  A. 
Regardless  of  the  type  of  estimation,  the  process  involved 
is  an  example  of  supervised  learning.   It  presupposes  a 
set  of  correctly  labeled  samples  from  which  the  estimates 
may  be  constructed.   Information  external  to  the  features 
being  measured  has  been  used  to  identify  the  members  of 
each  class. 

B.   DECISION  RULES 

The  previous  discussion  of  the  Bayesian  approach  to 
signal  classification  already  suggests  one  method  of  ob- 
taining a  decision  rule.   For  each  new  (unclassified) 
feature  vector  use  the  a-priori  probability  inferred  from 
the  taining  set  to  calculate  the  feature-conditional  pro- 
bability of  its  membership  in  each  class;  then  choose  that 
class  for  which  this  probability  is  greatest. 

Such  a  maximal  likelihood  criteria  is  intuitively 
satisfying  and  can  be  formulated  to  lead  to  classifications 
which  are  optimal  in  the  sense  of  minimizing  the  average 
risk  of  misclassif ication .   Unfortunately  a  naive  implemen- 
tation of  such  a  technique  leads  to  computational  inefficiencies 
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which  often  render  the  process  impractical  in  terms  of 
memory  and  time  constraints. 

More  frequently  the  decision  rule  chosen  in  a  practical 
pattern  recognition  problem  is  suboptimal  in  the  Bayesian 
sense,  providing  a  method  of  partitioning  the  feature-space 
of  the  samples  without  resorting  to  the  calculation  of  the 
maximal  likelihood  criteria.   These  suboptimal  schemes 
use  simplifying  assumptions  regarding  the  nature  of  the 
underlying  class-conditional  probability  distributions  to 
arrive  at  class  boundaries  which  are  relatively  simple  func- 
tions of  the  features  themselves.   If  there  is  little  overlap 
in  the  various  class-conditional  probability  density  func- 
tions, the  increase  in  the  expected  misclassif ication  rate 
as  a  result  of  the  suboptimal  decision  is  minimal. 

For  example,  consider  the  case  of  the  one-dimensional 
distributions  for  the  two  equiprobable  classes  shown  in 
Figures  la  and  lb.   Suppose  that  some  suboptimal  decision 
rule  results  in  decision  boundary 


B»  =  B   .  -  e  (II-4) 

opt 


The  shaded  areas  of  Figures  la  and  lb  indicate  the  resultant 
misclassif ication  rates. 

If  there  is  substantial  separation  of  the  class  distri- 
butions, the  differential  error  in  classification  as  a  result 
of  boundary  misplacement,  area  D,  is  negligible.   In  general, 
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.  (a)  Large  Overlap 


(b)  Small  Overlap 


FIGURE  1.   Misclassification  Rates  Due  ot  Suboptimal 
Decision  Boundaries 
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with  decreasing  overlap  of  unimodal  class  distributions, 
not  only  does  the  total  misclassif ication  error  decrease, 
but  so  does  its  sensitivity  to  the  placement  of  the  decision 
boundary.   As  a  practical  matter,  the  distinction  between 
a  five  percent  and  a  two  percent  misclassif ication  rate 
seems  unneccessarily  fine  if  the  probability  distributions 
were  originally  inferred  from  ten  samples  each. 

In  many  cases  the  sub-optimal  classification  schemes 
follow  rather  naturally  from  the  method  used  to  infer  the 
a-priori  distributions.   Unimodal  feature  distributions 
give  rise  to  classification  schemes  based  on  a  test  sample's 
minimal  "distance"  from  estimated  class  means.   More  complex 
feature  distributions  involve  the  use  of  a  distance  measure 
in  the  feature  space  to  establish  the  K  training  set  members 
closest  to  the  test  sample;  but  classification  based  on  the 
weighted  vote  of  these  K-nearest  neighbors,  seems  closer  in 
spirit  to  non-parametric  probability  estimation.   Both 
these  classification  methods  are  explored  more  fully  in 
Appendix  B. 
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III.   PRACTICAL  CONSIDERATIONS  IN  THE  IDENTIFICATION 
OF  NON-COOPERATIVE  SIGNALS 

The  foregoing  discussion  of  the  basic  techniques  of 
mathematical  pattern  recognition  is  intended  to  serve  as 
an  outline  of  approaches  to  the  general  problems  of  classi- 
fication.  Identification  of  the  originator  of  a  bauded 
signal,  while  bearing  considerable  similarity  to  Bayesian 
classification,  is  complicated  by  a  number. of  factors  not 
normally  arising  in  such  problems  as  the  identification  of 
handwritten  characters. 

In  many  of  the  classical  problems  of  mathematical  pat- 
tern recognition,  samples  are  obtained  from  a  statistical 
universe  which  is  well-defined  in  the  sense  that  the  number 
of  classes  to  which  the  samples  may  belong  is  finite  and 
previously  ascertained.   Additionally,  the  parameters  or 
features  used  to  typify  these  classes  are  normally  time- 
invariant.   The  set  used  for  training  is  accurately  classi- 
fied and  sufficiently  large  that  meaningful  statistical 
estimates  of  the  parameters  and  their  distributions  may 
be  obtained. 

Under  these  conditions,  the  processes  of   Section  I 
find  their  greatest  application.   Unfortunately,  they  are 
seldom  encountered  in  identifying  the  signals  of  a  non- 
cooperative  originator.   The  data  base  used  for  a  training 
set  does  not  normally  provide  arbitrarily  accurate  ground 
truth  classification.   In  practice,  even  the  number  of 
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actual  originator  classes  may  not  be  known.   In  many  cases, 
the  paucity  of  signal  intercept  precludes  more  than  the 
grossest  statistical  parameter  estimation.   The  parameters 
themselves  are  time-dependent  since  they  are  affected  by 
component  aging,  time-dependent  signal  propagation  condi- 
tions, and  unknown  external  maintenance  and  adjustment. 
The  advantages  of  long-term  averaging  to  obtain  statistical 
convergence,  i.e.  accurate  a-priori  estimates,  may  be 
obscured  by  short-term  perturbations  for  which  adequate 
correction  may  be  prohibitively  costly  or  undefinable. 

In  summary,  the  investigator  proposing  to  use  pattern 
recognition  techniques  to  identify  the  originator  of  hostile 
signals  is  faced  with  a  series  of  constraints  involving 
small  sample  size,  time-varying  signal  parameters  and  class 
probabilities,  inaccurate  ground  truth,  and  an  indeterminate 
number  of  classes  with  which  signals  may  be  associated. 
Such  constraints  impose  severe  limitations  upon  the  pattern 
recognition  techniques  which  may  be  successfully  applied  to 
the  problem.   The  following  sections  indicate  some  of  the 
procedures  which  deal  with  these  limitations. 

A.   SMALL  SAMPLE  SIZE 

The  difficulties  arising  from  small  sample  size  are  most 
immediately  apparent  when  one  attempts  to  increase  the  number 
of  features  in  order  to  facilitate  the  separation  of  classes. 
The  dominant  effect  is  one  of  inaccuracy  arising  from  under- 
sampling  the  required  d-dimensional  a-priori  probability 
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distributions.   When  this  undersampling  occurs,  the  most 
likely  result  is  that  increasing  the  number  of  features  be- 
comes counterproductive  to  "good"  classification.   In  the 
practical  multi-category  problem,  it  is  not  unusual  to 
encounter  an  initial  set  as  many  as  fifty  features  which 
the  designer  believes  to  be  useful  in  identification.   The 
number  of  taining  samples  is  often  too  small  to  allow 
meaningful  Parzen  window  or  other  non-parametric  density 
estimation  techniques  over  this  number  of  dimensions. 

An  alternative  procedure  of  characterizing  the  a-priori 
distribution  function  by  use  of  the  most  general  multi- 
variate normal  form  requires  estimation  of  a  "d  x  d"  non- 
singular  covariance  matrix.   This  imposes  an  algebraic  re- 
quirement of  at  least  d-M  samples.   However,  experience  in 
the  practical  estimation  of  covariance  matrices  suggests 
that  at  least  two  to  three  times  this  minimum  algebraic 
requirement  of  d+1  samples  is  preferable  [7,  11].   Frequent- 
ly the  number  of  available  training  samples  for  such  complete 
estimation  of  the  individual  class  covariance  matrices  is 
inadequate. 

Two  common  methods  of  correcting  this  difficulty  involve 
somewhat  arbitrary  assumptions  about  the  nature  of  the  class 
covariance  matrices.   One  possible  assumption  is  that  all 
classes  share  the  same  covariance  matrix  S,  which  is  then 
obtained  by  averaging  the  individual  class  covariance 
matrices  S . : 
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A  second  method  assumes  statistical  independence  of  the 
features  and  makes  all  off-diagonal  elements  zero,  regard- 
less of  evidence  to  the  contrary.   This  approach  preserves 
additional  detail  about  individual  class  structures  at  the 
expense  of  disregarding  highly  correlated  features.   It 
requires  only  the  calculation  of  a  single  feature  class 
means  and  variances,  as  indicated  in  Equations  A-l  and  A-2. 
Although  such  assumptions  are  almost  surely  incorrect,  they 
often  lead  to  better  classifier  performance  than  a  maximum 
likelihood  estimator  [7] . 

Small  sample  size  also  leads  to  complications  when  one 
attempts  to  estimate  the  error  rate  of  the  classification 
algorithm  finally  adopted.   If  one  chooses  to  estimate  the 
classifier's  error  rate  from  the  assumed  parametric  model, 
one  finds  the  result  is  optimistic  to  the  extent  that  the 
training  samples  are  peculiar  and  unrepresentative  [21] . 

An  empirical  approach  to  determining  error  rate  avoids 
this  problem  by  running  the  classifier  on  a  set  of  test 
samples.   Where  one  is  faced  with  a  small  number  of  available 
signals  for  which  ground  truth  is  accurate,  the  confidence 
which  may  be  placed  in  the  results  of  the  empirical  scheme 
is  marginal.   For  example,  if  two  errors  are  made  in  ten 
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test  samples,  one  may  predict  with  95%  confidence  that  the 
true  error  rate  of  the  classifier  lies  somewhere  between 
five  and  fifty-five  percent.   The  figure  of  ten  test  sam- 
ples is  not  unrealistic  if  one  uses  different  design  and 
test  samples  to  avoid  the  hazards  of  "testing  on  the  training 
data." 

B.   INDETERMINANCE  OF  CLASSES 

The  absence  of  a  predetermined  number  of  classes  presents 
additional  problems  which  must  be  addressed  in  any  practical 
identification  scheme.   The  first  problem  is  that  of  estab- 
lishing the  criteria  for  excluding  a  given  sample  as  a  member 
of  any  of  the  previously  determined  classes.   The  second 
problem  is  one  of  associating  or  "clustering"  the  resultant 
"unclassified"  signals  into  possible  new  classes. 

Within  the  formal  structure  of  minimal  risk  Bayesian 
classification,  one  may  introduce  the  concept  of  an  addi- 
tional class  which  corresponds  to  the  heuristic  category 
of  "don't  know."   More  frequently,  this  category  is  inter- 
preted as  that  portion  of  the  feature  space  for  which  no 
feature-conditional  class  probability  exceeds  a  given 
threshold.   That  is: 


p(C.|x)  <  X      for  all  i  (III-2) 


An  alternative  method  of  establishing  this  threshold 
for  multivariate  normal  distributions  is  to  specify  a 
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maximum  Mahalanobis  distance  from  the  established  class 
means . 


(i)  T   -1     (i) 
ri  "  (x-V  V  g±  <x-£   >  >  X  (HI-3) 


as  the  boundary  for  which  a  sample  will  be  included  in  a 
given  class. 

The  second  problem,  that  of  clustering  the  "unclassified" 
signals,  is  less  well  defined  from  the  standpoint  of  Bayesian 
analysis.   It  represents  one  of  a  group  of  problems  known 
as  unsupervised  learning.   Jarvis  and  Patrick  [12]  present 
several  techniques  by  which  such  clustering  may  be  performed 
and  illustrate  the  advantages  in  graphically  displaying 
the  clustering  process. 

One  such  clustering  technique  used  in  several  identifi- 
cation systems  is  illustrated  for  a  two-dimensional  case 
in  Figure  2.   Thresholds  based  on  Mahalanobis  distance  are 
established  for  each  cluster.   If  the  Mahalanobis  distance 
for  a  given  sample,  such  as  S..  ,  exceeds  the  outer  threshold 
distance  from  the  mean  of  any  known  cluster,  that  sample 
is  considered  to  represent  a  new  cluster  and  possibly  a  new 
signal  originator.   For  this  new  cluster,  arbitrary  dis- 
tance thresholds  based  on  "typical"  class  covariance  are 
assumed.   Subsequent  samples  (S2)  lying  within  the  inner 
threshold  are  used  to  update  the  class  statistics;  those 
samples  (S~)  lying  between  the  two  thresholds  are  tentatively 
associated  with  the  class,  but  do  not  update  the  class  statistics 
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FIGURE  2.   Guard  Zone  Clustering 
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It  is  during  the  clustering  procedure  that  the  man- 
machine  interface  becomes  a  critical  factor  in  the  pattern 
recognition  process.   An  equally  important  factor  in  the 
efficiency  of  the  clustering  algorithms  is  a  technique  for 
feature  selection  which  provides  an  adequate  basis  for 
representing  new  classes.   The  initial  investigator  has  the 
obligation  to  use  a  training  set  which  as  nearly  as 
possible  represents  the  entire  spectrum  of  possible  feature 
parameters . 

C.   INACCURATE  GROUND  TRUTH 

The  difficulty  of  obtaining  accurate  correspondence 
between  the  elements  of  the  initial  training  set  and  known 
signal  originators  represents  the  most  persistent  and  per- 
plexing dilemma  in  the  design  of  pattern  recognition  devices. 
Those  signals  for  which  there  is  a  high  degree  of  probability 
that  the  originator  can  be  determined  from  external  infor- 
mation may  not  adequately  represent  the  entire  spectrum 
of  signals.   Conversely,  the  originators  of  a  more 
representative  group  of  signals  may  not  be  known. 

One  may  attempt  to  choose  self  consistent  classes 
through  the  use  of  a  clustering  algorithm  and  the  applica- 
tion of  probability  estimator  to  information  confirming 
class  membership,  a  technique  best  described  as  "learning 
with  a  probabilistic  teacher."   This  approach,  discussed 
more  fully  by  Cooper  [4] ,  is  still  rather  exploratory  in 
nature  and  is  implemented  at  considerable  cost  in  timeliness 
and  computational  efficiency. 


27 


Since  the  empirical  methods  by  which  any  new  classifi- 
cation performance  is  evaluated  must  use  as  ground  truth 
the  identifications  of  an  older  system,  indications  of  any 
improved  performance  in  the  new  system  are  perforce  intui- 
tive.  This  is  the  crux  of  the  circular  dilemma  which  arises 
when  one  is  forced  to  estimate  not  probability  of  error,  but 
probability  of  disagreement.   In  the  event  of  disagreement, 
which  system  is  right? 
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IV.   FEATURE  SELECTION  IN  WAVEFORM  ANALYSIS 

One  important  question  regarding  the  process  of  waveform 
identification  has  not  yet  been  addressed.   All  of  the  tech- 
niques described  so  far  have  indicated  that  one  commences 
this  process  with  a  particular  set  of  feature  measurements. 
In  the  case  of  a  waveform  z(t),  a  sufficiently  general 
feature  vector  might  be  the  sequential  samples  of  the 
continuous  signal,  i.e.  {z (t, ) , z (t«) . . . ,z (t  )}.   If  these 
samples  occur  at  greater  than  the  Nyquist  rate,  one  may 
reasonably  assume  that  they  constitute  an  adequate  descrip- 
tion of  the  signal.   One  could  then  proceed  directly  to 
the  calculation  of  n-dimensional  probability  distributions 
for  each  class,  blithely  oblivious  to  the  fact  that  n 
may  be  on  the  order  of  500  or  more. 

In  light  of  the  advantages  in  the  choice  of  a  small  number 
of  mutually-independent  well-separated  features  for  the 
process  of  pattern  identification,  it  is  normally  advisable 
to  seek  some  mapping,  F,  of  the  signal  space  Z  into  a 
pattern  space  X.   That  is,  we  wish  to  find  a  preprocessing 
transformation  so  that  the  pattern  space  X  =  F(Z)  satisfies 
the  following  objectives: 

(1)  low  dimensionality 

(2)  retention  of  sufficient  information 

(3)  enhancement  of  distance  in  pattern  space  as  a  measure 
of  the  similarity  of  physical  patterns 
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(4)   comparability  of  features  among  samples. 
Note  that  1  and  2  imply  elimination  of  redundant  information 

In  addition  to  the  above  objectives  one  might  also  wish 
to  obtain  a  pattern  space  representation  adequate  for  the 
construction  of  new  classes  in  an  unsupervised  learning 
mode.   To  this  end,  it  is  convenient  that  the  individual 
components  of  the  pattern  space  have  some  natural  interpre- 
tation which  might  provide  qualitative  information  about 
the  underlying  causes  of  class  difference. 

As  a  first  step  toward  eliminating  redundant  informa- 
tion in  a  pattern  space  satisfying  the  above  criteria,  it 
is  frequently  desirable  to  represent  the  sampled  waveform 
as  a  linear  combination  of  orthonormal  functions .   Two 
commonly  used  orthonormal  expansions  satisfying  slightly 
different  training  set  mean  square  error  criteria  are  dis- 
cussed in  Appendix  C.   The  coefficients  of  these  functions 
may  then  be  used  as  features  of  the  signal.   The  choice  of 
orthonormal  functions  used  in  the  expansion  however,  need 
not  be  one  which  necessarily  minimizes  a  form  of  mean 
square  error  in  the  representation  of  the  signal.   If  one 
is  willing  to  accept  a  slightly  greater  amount  of  error  in 
modeling  the  training  set  waveforms,  then  one  possible 
technique  is  to  use  a  set  of  orthonormal  functions  for 
which  it  is  computationally  efficient,  as  it  is  with  the 
fast  Fourier  transform,  to  obtain  coefficients.   By  trun- 
cating the  number  of  coefficients  one  may  achieve  some 
reduction  of  dimensionality. 
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The  features  derived  from  the  ortho-normal  projections 
of  the  sampled  signal  described  in  Appendix  C  represent 
linear  transformations  of  the  sampled  waveform  data.   For 
the  purposes  of  signal  identification,  this  representation 
may  not  be  adequate. 

Since  one  of  the  goals  in  this  application  of  pattern 
recognition  techniques  is  to  assist  a  human  analyst  in 
making  an  identification,  it  is  particularly  desirable  that 
the  resultant  features  derived  from  preprocessing  have  some 
degree  of  natural  interpretation.   The  optimal  linear  trans- 
formations discussed  do  not  always  lend  themselves  to  easy 
extension  in  an  unsupervised  learning  mode. 

Features  which  are  more  directly  related  to  an  assumed 
model  of  the  signal  generating  process  and  which  are  often 
easily  implemented  as  a  measurement  procedure,  may  in  fact 
involve  non-linear  transformations  of  the  waveform  data.   In 
selecting  features  of  this  type,  the  prior  experience  of 
the  data  analyst  and  the  details  of  the  signal  generating 
model  serve  as  a  focus  for  the  development  of  non-linear 
preprocessing. 

Quite  commonly,  the  initial  feature  space  of  the  signal 
will  include  components  derived  from  readily  interpretable 
data  transformations,  both  linear  and  non-linear  in  nature. 
This  feature  set  is  normally  too  large  for  convenient  compu- 
tation and  some  of  the  features  thus  obtained  may  be  depen- 
dent.  The  techniques  of  linear  dimensionality  reduction 
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described  in  Appendix  C  may  be  applied  to  this  set  of 
features,  particularly  when  graphical  display  of  the  indi- 
vidual signal  is  required.   In  order  to  provide  a  feature 
set  with  dimensionality  low  enough  to  insure  computationally 
convenient  classification  and.  general  enough  to  allow  the 
establishment  of  new  class  clusters,  it  is  both  sufficient 
and  desirable  to  reduce  dimensionality  while  retaining  the 
identity  of  individual  features  by  selecting  those  n  of 
m  features  which  lead  to  the  optimal  separation  of  classes. 
The  selection  process  then  requires  both  the  adoption  of  a 
criterion  by  which  class-separability  may  be  estimated  and 
the  development  of  a  feature  search  algorithm.   Both  of 
these  topics  are  discussed  in  Appendix  C. 
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V.   PARAMETER  MAPPING  OF  RASTER  SCAN  DISPLAYS 

As  a  practical  application  of  the  signal  classification 
techniques  previously  described,  the  author  conducted  a 
series  of  investigations  which  used  as  a  data  base  the  raster 
scan  displays  of  bauded,  frequency-shift-keyed  (FSK)  sig- 
nals transmitted  by  several  originators.   The  signal  data 
base,  many  of  the  associated  processing  algorithms  and  the 
processing  equipment  itself  were  made  available  through  the 
cooperation  of  Electromagnetic  Systems  Laboratories  in 
Sunnyvale  California.   The  data  base  and  supporting  facili- 
ties are  the  results  of  the  company's  development  of  the 
Parameter  Encoder  in  a  Navy-sponsored  program  for  signal 
measurement  and  identification. 

The  objective  of  the  author's  investigation  was  to  iso- 
late from  the  raster  scan  data  a  small  set  of  features  which 
might  prove  useful  in  automated  signal  identification.   As 
a  part  of  the  existing  clustering  and  identification  tech- 
nique, the  raster  scan  pattern  of  a  given  signal  was  pre- 
sented to  the  system  operator  for  visual  comparison  with 
patterns  of  previously  identified  signals.   It  was  claimed 
that  such  comparison  was  of  value  in  the  final  stages  of 
the  identification  process,  when  the  classification  of  the 
signal  in  question  had  been  narrowed  to  a  few  possibilities. 
Display  of  the  information  necessary  to  produce  the  required 
number  of  raster  scan  patterns,  however,  consumes  a  significant 
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amount  of  computer  time.   The  interests  of  improved  effi- 
ciency, coupled  with  the  author's  desire  to  quantify  a 
process  open  to  considerable  latitude  in  operator  judgment/ 
provided  the  motivation  for  the  project. 

In  obtaining  a  reasonable  initial  set  of  piirameters  for 
consideration  by  this  thesis,  a  study  of  the  signal  generating 
process  and  the  characteristics  of  the  raster  scan  display 
led  to  a  slight  modification  of  the  normal  minimal  phase 
raster  representation.   The  resulting  "phase-unwrapped" 
raster  display  was  the  basis  of  a  series  of  measurements 
related  to  clockrate,  transient  phenomena,  and  possible 
data- dependent  or  unintentional  external  rate  modulation. 
A  model  of  the  signal  process,  the  raster  scan  display,  and 
the  methods  used  to  measure  the  initial  set  of  features 
are  discussed  in  the  following  sections. 

A.   SIGNAL  GENERATION  PROCESS 

Since  the  raster  scan  display  is  based  essentially  on 
zero  crossing  information,  it  represents  a  non-linear  trans- 
formation of  the  incident  signal.   It  is  convenient  to  refer 
to  a  model  of  the  signal  generating  process  in  an  attempt 
to  justify  use  of  such  a  transformation.   From  this  model 
the  qualitative  effects  of  the  various  steps  in  the  genera- 
tion of  the  raster  scan  pattern  can  be  estimated.   Fig.  (3) 
provides  an  outline  of  the  process  which  leads  to  the  raster 
scan  display.   All  of  the  signal  data  base  information  avail- 
able for  the  feature  measurements  proposed  by  this  thesis  is 
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contained  in  the  transition  time  measurements  provided  by 
the  edge  detector. 

Of  all  the  qualitative  effects  upon  transition  time 
measurements,  those  associated  with  the  originator's  clock 
and  bit-stream  generator  are  potentially  of  greatest  value 
to  a  raster-scan-based  identification  scheme.   Class-to- 
class  differences  in  the  basic  clock  rate  are  certainly  of 
interest,  as  are  clock  rate  variations  due  to  component 
instability  and  unintentional  external  rate  modulation. 
If  the  originator's  clock  is  not  continuously  running, 
transient  phenomena  associated  with  clock  turn-on  may  be 
apparent  in  the  raster  display.   The  binary  bit  stream 
generator  (data  modulator)  and  non-linearities  in  the  FSK 
modulator  itself  may  produce  an  unequal  mark-space  duty 
cycle  or  data  dependent  mark-space  asymmetries  which  appear 
as  bias  in  the  raster  display.   Such  effects  should  produce 
features  of  the  raster  display  which  serve  to  characterize 
the  originator.   To  isolate  them  effectively,  however,  may 
require  the  use  of  non-linear  waveform  data  transformations. 

Other  parts  of  the  signal  generating  process  may  intro- 
duce qualitative  effects  which  tend  to  obscure  those  produced 
by  the  source.   For  example,  noise  and  signal  attenuation 
in  the  propagation  path  may  introduce  FM  spikes  and  anoma- 
lous transition  times.   Non-linearities  and  zero-level  off- 
set in  the  FM  discriminator  may  produce  additional  bias  in 
the  measured  mark  and  space  durations.   The  waveform  sampling 
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procedure  itself  introduces  an  unavoidable  time  quantization 
error  which  may  or  may  not  be  uniformly  distributed  on  a 
short  term  basis.   If  this  error  is  uniformly  distributed, 
one  can  expect  a  transition  time  error  variance  a„  where 


a  =  ££  V-] 

°T    12  v  " 


and  where  At  is  the  sampling  period. 

B.   RASTER  SCAN  DISPLAY 

The  use  of  a  raster  scan  display  arises  rather  naturally 
in  the  attempt  to  represent  the  fine  details  of  clock  rate 
variation.   The  axes  of  the  display  correspond  to  coarse  and 
fine  divisions  of  time  in  the  following  manner.   Suppose 
some  event,  such  as  a  mark-space  transition,  is  assumed  to 
occur  only  at  times  characterized  by  some  integer  multiple 
of  a  nominal  period  T.   The  actual  time  of  occurence,  t,  may 
then  be  represented  by 

t  =  nT  +  z 

-T/2  <  z  <_  T/2;   n  integer  V-2 

By  equating  the  horizontal  axis  of  the  raster  display  to 
the  actual  time  variable,  t,  and  the  vertical  axis  to  the 
fine  time  z,  any  arbitrary  time  may  be  represented  as  a 
point  on  the  sawtooth  waveform  indicated  in  Figure  4b.   The 
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(a)  Phase  and  Phase  Perturbation  vs.  t 


(b)  Raster  of  Events  £,-£, 


FIGURE  4.  Relation  of  Phase  Perturbation  to  Raster  Display 


raster  display  has  the  effect  of  mapping  the  points  t. 

in  a  one-dimensional  region  into  the  points  (z(t.),t.) 

of  a  two-dimensional  region.   Since  the  fine  time  variable, 

z,  represents  a  lag  or  lead  relative  to  the  nominal  period, 

one  can  relate  the  points  of  the  raster  scan  display  to 

the  samples  of  some  continuous  phase  perturbation. 

For  example,  consider  the  continuous  phase  function 
illustrated  in  Figure  4a.   The  function  is  comprised  of  a 
linear  part  and  a  perturbation 


$(t)  =  2fr(i  t  +  z(t))  V-3 


where  T  is  the  nominal  period  of  the  raster  scan  display 
at  times  t...   If  samples  of  this  function  are  taken  when 


<&(t  -  )  =  2nir    n  =  0,1,2.  .  .  V-4 

n 


and  if 

|z(t) |  <  T/2   for  all  T  V-5 

then  the  minimal  time  difference  from  the  expected  sample 
time  is  proportional  to  the  phase  perturbation.   The  raster 
scan  points  of  Figure  4b  represent  non-uniform  samples  of 
this  continuous  perturbation.   In  this  sense  one  can  repre- 
sent the  events  occuring  at  times  (t, ,  t^  ...)  as  samples 
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(z(t..)  z(t2)  ...)  of  a  continuous  phase  proportional  function 
z(t)  where 


z (t  )  =  t   -  n    T 
n     n    max 


where 


"max  "  ^   +  ¥  V~6 

t.       , 

i.e.  the  gjreatest  integer  less  than  [-=-  +  y]  . 

Figures  5a  through  5d  illustrate  the  raster  patterns 
of  events  t.  of  arbitrary  starting  phase  and  periodicity 
different  from  the  nominal  period  of  the  raster  scan.   If 
the  nominal  raster  period  coincides  with  the  actual  periodic- 
ity, as  in  Figures  5a  and  5b,  the  slope  of  the  line  joining 
points  of  the  raster  pattern  is  zero.   Similarly,  if  the 
periodicity  of  the  series  of  events  differs  slightly  from 
the  nominal  raster  period  by  a  small  amount  e,    the  raster 
pattern  may  be  characterized  by  a  line  whose  slope  is  given 
by: 


z9  -  zn 

_± L  -   — £ —  V-7 

t2  -  tl    T  +  £ 


as  shown  in  Figures  5c  and  5d. 

In  the  signal  display  of  the  Parameter  Encoder  both 
mark-space  and  space-mark  transition  times  are  plotted  on  a 
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raster  scan  whose  nominal  period  is  approximately  one-half 
the  time  for  a  single  mark-space  cycle,  that  is  raster 
period  is  equal  to  the  baud  length.   The  alternate  points 
of  the  raster  scan  exhibit  a  bias  term  if  the  mark  and 
space  durations  are  not  equal.   If  the  period  of  the  entire 
mark-space  cycle  remains  constant,  lines  connecting  alter- 
nate points  of  the  raster  scan  have  the  same  slope,  but  they 
are  vertically  offset,  by  the  difference  in  mark  and  space 
durations.   Figure  6  illustrates  this  condition. 

The  effect  of  the  uniform  sampling  of  the  waveform  and 
the  subsequent  time  quantization  of  events  is  apparent  when 
the  nominal  period  of  the  raster-  scan  is  on  the  order  of 
ten  times  that  of  the  sampling  period.   Figures  7  and  8 
illustrate  the  effect  of  quantization  noise  in  two  raster 
scan  patterns  of  a  slowly  varying  signal,  one  of  which 
uses  a  nominal  period  that  is  is  an  integer  multiple  of 
the  sampling  time.   This  quantization  phenomena  was  evident 
in  the  raster  displays  of  the  signal  data  base  used  in  this 
thesis  since  the  nominal  baud  length  is .on'  the  order  of 
twenty  times  the  sampling  period. 

The  abrupt  transition  in  Figures  7  and  8  also  serves  to 
illustrate  the  "aliasing"  effect  which  may  occur  in  a  minimal 
phase  representation  if  the  time  difference  between  the 
actual  and  the  nominal  event  exceeds  one-half  the  nominal 
period.   When  this  aliasing  occurs  the  points  plotted  on 
the  raster  scan  are  no  longer  proportional  to  the  assumed 


42 


CO 

< 


o 

< 
CQ 


.J 
O 

o 


o 

> 

(0 

c 

•H 
CO 

13 

cu 
to 

fCS 
•H 
« 


rrj 


I  I 

>* 


I     co 

< 

0Q 


I" 


I 


X 


CO 

tf 

•H 

PQ 

iH 

fd 

£ 

Ln 

•H 

CO 

<a 

rH 

0) 

KJ 

*a 

c 

p 

CT> 

03 

•H 

CQ 

CO 

U-! 

1)' 

O 

(U 

CO 

>i 

m 

a( 

•H 

rH 

cq 

0* 

CO 

in 

■H 

0 

Q 

^ 

U 

0) 

0) 

4J 

-J-' 

to 

CO 

cci 

td 

« 

& 

XI 


w 

D 
O 
H 


H 


.^ 


I 


43 


*. 


-p 

<H 

t= 

fM 

G 

'H 

10 

< 

u 

£ 

to 

0 

-p 

•H 

■H 

-P 

£ 

N 

rrj 

3 

K 

♦Q 

H 

o 

in 

3 

CM 

• 

-P 

H 

II 

1! 

0) 

fr 

< 

IH 

ro 


(M 


wo='«      - ~ 


o 


>1 

rrj 

iH 

o. 

w 

•H 

P  i- 

d 

5-1 

> 

(D 

5h 

4J 

o 

CO 

-p 

fd 

c 

« 

H 

d 

c 

•H 

C) 

•rH 

CO 

-P 

•P 

rd 

o 

N 

CD 

-H 

HH 

P 

MH 

G 

W 

fd 

2 

tn 

O 

G 

•h 

X 

w 

rd  o 

•H 

• 

rH 

CO 

rfj 

rH 

13 

II 

C 

fd 

-d 

0 

C 

•H 

0 

5-1 

■H 

0) 

+J 

CLl 

rd 

N 

rH 

-H 

rd 

P 

C 

c 

•H 

fd 

e 

3 

0 

OS 


H 
OS 
D 
O 

M 
[=4 


^ 


44 


-p 

i= 
CM 


C 

•H 
CO 

< 

II 

C 
O 

•H 
•P 

fd 
,Q 

U 

-P 
H 


to 
-P 
•H 

o 


< 


N 

a 

IT) 


m 


ro 


<M 


>1 

rd 

rH 

a 

CO 

•H 

Q 

H 

0) 

-P 

CO 

fd 

to 

tf 

-p 

•H 

C 

C 

•H 

E3 

CO 

C 

•P 

0 

U 

•H 

<D 

■P 

in 

fd 

HA 

N 

w 

•H 

■P 

tn  c 

C 

fd 

•H 

3 

CO 

a 

fd 

•H 

rH 

o 

< 

CM 

•o 

II 

c 

fd  HJ 

0 

G 

•H 

o 

IM 

•H 

a) 

■P 

PU 

fd 

N 

rH 

•H 

fd 

•P 

r, 

c 

•H 

fd 

e 

3 

0 

a^ 


! 
-4- 

o 


r* 


CO 


O 
M 
P4 


45 


continuous  perturbation.   In  order  to  maintain  proportionality 
it  is  preferable  to  represent  the  time  of  events  by  a  minimal 
phase  different  plot.   Such  a  procedure  performs  a  certain 
amount  of  "phase  unwrapping"  by  assuming  that  the  absolute 
raster  time  difference  between  two  successive  points  never 
should  exceed  T/2 .   Figure  9  illustrates  the  effect  of  this 
minimal  phase  difference  algorithm.   Since  each  segment  of 
the  sawtooth  pattern  raster  scan  may  be  extended  in  either 
forward  or  backward  in  time  the  minimal  phase  difference 
representation  uses  this  extension  to  represent  those  points 
for  which  the  minimal  phase  difference  is  less  than   T/2  . 
In  Figure  9,  those  points  occurring  in  the  heavily  lined 
region  indicate  the  normal  raster  scan  pattern  while  a 
continuous  line  joins  the  points  of  the  minimal  phase 
difference  pattern.   Since  some,  of  the  "aliasing"  associated 
difficulties  in  the  measurement  of  features  may  be  eliminated 
by  use  of  a  minimal  phase  difference,  this  representation 
was  included  as  the  initial  step  in  the  measurement  process. 
Figure  2  6  outlines  the  flow  of  the  algorithm  actually  used 
in  the  measurement  program. 

C.   RASTER  SCAN  DESCRIPTIVE  PARAMETERS 

The  first  efforts  of  the  author's  investigation  of  the 
raster  scan  data  were  directed  towards  obtaining  a  rather 
large  set  of  descriptive  parameters,  particularly  those 
which  seemed  useful  in  the  quantitative  representation  of 
the  effects  anticipated  from  a  study  of  the  model  signal 
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generating  process.   Manipulations  of  the  transition  time 
data  to  provide  information  about  the  nominal  clock  rate, 
it's  transient  variations  and  external  modulations,  as  well 
as  source-associated  mark-space  bias  and  data-dependent 
transition  time  anomalies  were  the  areas  of  particular 
interest.   The  parameters  considered  included  representa- 
tions of  all  or  part  of  the  signal's  raster  phase  by  poly- 
nomial mean  square  fits,  Laguerre  polynomials,  and  Fourier 
power  spectrum  terms,  as  well  as  measures  of  mean  bias, 
intrinsic  and  total  signal  variance.   The  techniques 
employed  to  calculate  these  parameters  are  described  in 
the  f ollowing  sections . 

-*•  •   Mean  Bias^and  the  Effect  of  Transients 

Although  several  sources  of  apparent  mark-space  bias 
other  than  those  associated  with  the  originator  have  been 
indicated  in  the  study  of  the  process  model,  the  measure- 
ment of  this  parameter  may  provide  an  adequate  estimate  of 
the  source-associated  component.   More  recent  studies  of 
new  signal  data  for  which  the  bias  effects  not  associated 
with  the  originator  are  believed  to  be  negligible,  indicate 
that  such  source-dependent  bias  not  only  exists  but  can  in 
fact  be  time-varying. 

The  simplest  estimate  of  mean  bias  may  be  obtained 
by  the  straightforward  averaging  of  the  raster  phase  differ- 
ence of  alternate  points  on  the  raster  scan.   For  the  raster 
scan  display  of  the  alternate  positive  and  negative  transitions, 
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-        +        -  +  A 

z  (t.),  z  (t~) ,  z  (t^)  .  ..,  z  (t?  ),  the  mean  bias  b  is 

given  by: 


1   n    + 
b  =  i  E   z  (t,4)  -  z"(t„.  ,) 
n  .  ,      21       zi-i 
1=1 


1  n   +  i  n   - 

~     I      z  (t9.)  -  i  E  z  (t     ) 

n  .  ,      21  n  .  .     2i-l 

i=l  i=l 


=  z   -  z  V-8 


which  is  the  difference  in  the  mean  raster  phase  of  the  two 
bias  conditions. 

One  difficulty  which  arises  from  this  representation 
of  mean  bias  occurs  as  a  result  of  clock  transient  behavior. 
If  at  any  point  t.,  raster  phase  z(t.)  differs  from  the 
raster  phase  of  the  starting  point  by    an  amount  approxi- 
mately equal  to  one-half  the  nominal  raster  period,  the 
cyclic  representation  (minimal  phase)  of  the  raster  display 
introduces  an  "aliasing"  error  into  the  mean  bias  estimate 
of  Equation  V-B. 

To  eliminate  such  errors  in  the  bias  estimate,  one 
must  use  the  minimal  phase  difference  representation  discussed 
previously  for  each  bias  condition. 

If  the  transient  behavior  of  the  clock  rate  is  suffi- 
ciently violent,  even  the  minimal  phase  representation  may 
not  be  adequate.   For  example,  assume  that  the  instantaneous 
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originator  clock  frequency  f,  is  characterized  by  turn-on 
from  a  dead  stop: 


f  =  f   [1  -  exp(-at)  ]  V-9 

s  ^> 


where  a  is  the  reciprocal  of  turn-on  constant,  t.   The 
instantaneous  phase  <f>(t)  in  cycles  is  obtained  by  integrating 
f  from  0  to  t: 


f 

<Mt)  =  (f  c  t  +  -M  [-1  +  exp  (-at)]}    V-10 
So      a 


If  transitions  occur  when  <J>(t)  =  0,1,2,...,  (Figure  10a), 

the  raster  scan  points  of  Figure  10b  are  obtained.   It  should 

be  noted  that  the  transition  times  are  non-uniform  samples 

of  the  actual  phase.   Similar  computations  may  be  performed 

for  the  raster  phase  of  a  clock  whose  transient  frequency 

behavior  is  characterized  by  a  shift  from  a  free  running 

frequency,  f  ,  to  the  steady  state  frequency,  f 

o  s  s 


f   =  f   exp  (-at)  +  f  _  [1  -  exp  (-at)] 
o     o  s  s 


■  [fo  -  fss]  e*P  (~at)  +  fss  V"n 


[f  -  f  ] 

<j)(t)  =  — £— ^~  [exp  (-at)  -1]  +  f    t     V-12 


The  above  results  suggest  that  the  raster-phase 
pattern  or  a  minimal  phase  difference  pattern  may  be 
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extremely  misleading  during  the  transient  period.   It  is 
evident  that  the  predicted  raster  phase  (in  cycles)  may 
differ  from  the  actual  steady  state  phase  by  an  amount 
equal  to 


f   -  f      T    -  T 
o ss  _   ss o 

a         T   T    T 
o   ss 


which  may  be  considerably  greater  than  one  cycle.   For 

example,  if  T  =  ,9t   and  t  =  20  T   ,  then  difference  in 
r    ■ '  o      ss  ss' 

steady  state  phase  may  be  more  than  two  complete  cycles. 

Using  the  same  time  constant  x,  turn-on  from  a  dead  stop 

(f   =  0)  results  in  a  steady  state  phase  difference  of 

twenty  cycles  (from  nominal) .   The  raster  scan  display  of 

Figure  11  serves  to  illustrate  the  transient  turn-on  effects 

observed  when  t  =  10  T   .   In  such  a  case  the  anti-aliasing 

ss  ^ 

capability  of  the  phase  unwrapping  process  is  severely 
strained  and  it  would  be  wise  to  avoid  measurement  of  mean 
bias  in  this  transient  region. 

2 .   Pattern  Variance  and  Intrinsic  Variance 

Pattern  variance  and  intrinsic  variance  measure 
the  extent  to  which  the  raster  pattern  is  explained  by  the 
nominal  clock  rate  and  the  minimal  phase  difference  represen- 
tation.  Large  values  of  total  pattern  variance  indicate 
that  the  nominal  raster  chosen  does  not  completely  explain 
the  pattern.   Large  values  of  the  intrinsic  variance  suggest 
that  a  minimal  phase  difference  representation  is  not 
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adequate,  perhaps  as  a  result  of  violent  transients  or 
noise. 

These  two  measures  are  averaged  over  points  of  both 

(1)       (2) 
bias  conditions  z    and  z    corresponding  to  the  alternate 

positive  and  negative  going  transitions.   The  total  variance 

V™  is  obtained  from  the  following  calculation: 

9  N.  N. 

VT   "  .\   [N^T  ,E,  Zi      +  N.(N.-1)(.Z.,  Zj   >  3     V"14 
1=1    l    3=1   J         li     3=1   J 


Intrinsic  variance  (vT)  is  estimated  by: 

?  Ni/2 

vT2  =   Z  -L  Z    (z(l)  _  2(i)  ;2 
1     i=l  Ni  j=l    ^     ^l' 

2    3   N±  (z(i)~  7^i2  V-15 

1=1    l  3=2 

where  z.     refers  to  the  jth  point  of  the  set  of  either 
3 

positive  or  negative  going  transitions. 

These  two  terms  may  be  particularly  useful  as  a 
measure  of  signal  degradation.   They  may  be  used  as  a  measure 
of  non-parametric  correlation  since  the  fraction  of  total 
variance  "explained"  by  the  intrinsic  variance  is  given  by: 

2 

A(z)  =1 ~  V-16 

VT 
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One  estimate  of  non-parametric  correlation  is: 

B(z)  -   VMz)   •  V-17 

Small  values  of  B(z)  suggest  that  the  apparent  behavior  of 
the  raster  pattern  is  not  well  characterized  by  a  relatively 
smooth  function  of  time. 

3 .   Polynomial  Mean  Square  Fits 

In  the  hope  of  quantifying  raster  pattern  behavior 
which  system  operators  had  described  as  "curve  up,"  "curve 
down,"  flat,"  or  complex,"  a  fourth  order  polynomial  mean 
square  fit  to  the  minimal  phase  difference  raster  data  over 
the  entire  signal  period  was  attempted  for  each  bias  condi- 
tion.  The  fourth  order  fit  was  initially  believed  to  be 
capable  of  representing  both  transient  and  steady  state 
phenomena.   Unfortunately,  even  the  minimal  phase  difference 
representation  of  the  data  could  not  eliminate  aliasing 
problems  in  those  displays  where  apparently  large  initial 
transients  occurred.   As  an  alternative,  a  second  order  fit 
over  that  data  lying  outside  the  region  of  worst  transient 
behavior  was  eventually  incorporated  into  the  measurement 
program. 

The  coefficients  of  the  second  order  fit  also  serve 
the  purpose  of  reducing  the  sensitivity  of  the  Fourier  trans- 
form data  to  phase  offset  and  endpoint  discontinuities  which 
result  from  minor  adjustments  in  the  nominal  raster  period 
of  the  display. 
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As  an  example,  consider  the  periodic  Fourier  series 
transform  of  a  pattern  characterized  by  a  line  with  non-zero 
slope.   As  a  result  of  endpoint  discontinuity  in  a  periodic 
representation,  the  odd  Fourier  series  components  exhibit 
amplitudes  proportional  to  the  magnitude  of  the  discontinuity 
and  decreasing  in  frequency  as  — . 

It  has  previously  been  established  that  transient 
phenomena,  zero  point  offset,  and  differences  between  the 
nominal  and  actual  steady  state  period  can  produce  patterns 
with  both  offset  and  non-zero- slopes .   If  the  raster  pattern 
shows  significant  curvature,  it  is  difficult  to  justify  any 
particular  choice  for  the  nominal  period.   Any  choice  of 
nominal  period  creates  its  own  contributions  to  the  Fourier 
series  components.   Since  the  nominal  raster  period  of  the 
display  is  generally  the  result  of  an  automatic  measurement 
procedure,  and  since  it  represents  one  of  the  signal  features 
currently  used  by  the  system's  identification  algorithm,  the 
interests  of  comparability  were  best  served  by  partially 
correcting  for  these  discontinuities  in  terms  of  the  mean 
square  fit  coefficients.   Only  the  difference  between  the 
raster  pattern  and  the  polynomial  mean  square  fit  to  each 
bias  condition  was  used  as  input  to  the  fast  Fourier  transform 
calculation. 

4 .   Fourier  Power  Spectrum 

In  order  to  observe  the  effects  of  possible  uninten- 
tional external  modulation  and  periodic  variation  of  the  clock 
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rate,  the  use  of  Fourier  power  spectrum  components  was  pro- 
posed.  The  calculation  of  these  terms  presented  an  intriguing 
problem,  since  the  minimal  phase  raster  data  represent  non- 
uniform time  samples  of  a  presumably  continuous  signal.   The 
non-uniformity  of  sampling  time  arises  from  both  the  clock 
rate  variation  and  the  fact  that  the  transition  times  them- 
selves and  data  modulated.   The  presence  of  one  transition 
per  baud  is  not  the  case;  transitions  occur  at  various 
integer  numbers  of  originator  clock  cycles. 

The  problem  of  non-uniform  sampling  is  not  too  seri- 
ous if  one  intends  to  obtain  only  a  few  Fourier  power  series 
components.   The  techniques  for  estimating  amplitude  phase 
and  frequency  of  the  principal  Fourier  components  has  been  a 
topic  in  recent  literature  [22].   One  can  easily  construct 
Fourier  series  of  terms  which  are  orthogonal  on  the  interval 
[0,L]  by  considering  only  the  frequency  terms  which  are 
periodic  in  L. 

.27m, 

°=         l-y-— t 

z(t)  =   Z   c  e   Jj  V-18 

n»—  n 

These  complex  components 

.  2Tm  . 
L      ~1_f —  t 
c_  =  /  z(t)e   L  dt  V-19 

n    0 

may  be  approximated  step-wise 
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-i|22_t. 

n  t     1 

max  n 

c   =   S   z(t.)  e    max   (t,-t.  .)  V-20 

n    i=2     i  i   i-l 


or  by  numerical  integration  involving  higher  order  interpo- 
lation.  Although  numerical  integration  of  the  above  type  is 
possible,  all  frequency  components  up  to  the  nominal  clock 
frequency  were  potentially  of  interest  as  features.   An 
alternate  form  of  calculating  the  magnitudes  of  the  power 
spectral  components,  using  the  fast  Fourier  transform  was 
adopted  for  computer  implementation . 

For  each  bias  condition,  this  procedure  typically 
used  about  75  non-uniformly  spaced  raster  data  points  obtained 
outside  the  transient  region.   For  each  condition,  12  8 
sample  points  uniformly  spaced  over  the  same  time  interval 
were  chosen.   The  raster  phase  at  these  points  was  obtained 
by  an  interpolation  procedure  which  performed  a  cubic  poly- 
nomial fit  to  four  minimal  phase  difference  raster  points, 
two  on  either  side  of  the  interpolation  point.   Since  the 
fast  Fourier  transform  routine  used  accepts  complex  values 
as  input,  the  interpolated  functions  for  positive  going  and 
negative  going  transitions  were  stored  as  the  respective 
real  and  imaginary  parts  of  this  input.   The  complex  input 
function  may  be  expressed  as: 


y  =  f  +  ig  V-21 


where  the  nth  Fourier  component  of  the  real  functions  f  and 
g  are: 

Fn(f)  -  ax  +  i  bx 

Fn(g)  =  a2  +  i  b2  V-22 

The  linearity  of  the  Fourier  transform  and  the  requirement 
that  for  a  real  function 


F-n(f)  =  Fn(f)  V"23 


implies  that  the  nth  Fourier  transform  components  of  both 
functions  f  and  g  may  be  recovered  from  (F  (y)  and  F_  (y) ) , 
the  real  and  imaginary  Fourier  components  of  the  complex 
input  [10]  .   Additionally  it  can  be  shown  that: 


[Re[Fn(y)]]2  +  [Im[Fn(y)j]2  +  Re (F_n (y) ) ] 2  +  [Im[F_n (y) J ] 2 
=  2(aq2  +  b  2  +   a22  +  b22)  V-24 


The  Fourier  transform  information  used  as  features  in  this 
thesis  was  the  squared  magnitude  of  the  frequency  components 
summed  over  both  bias  conditions  in  the  manner  shown  above. 

The  Figures  12-22  show  the  behavior  of  the  time 
quantization  interval,  the  minimal  phase  difference  interpola- 
tion, the  second  order  polynomial  fit,  and  the  Fourier  transform 
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features  on  a  signal  whose  possible  transition  times  t 

n 

are  given  by 


t   =  (nT  +  Asin  27rfnT  +  B  exp(-anT)) 
n  c 


The  data  points  were  selected  from  these  possible  transition 
times  by  a  simulated  data  modulation  which  repeated  itself 
about  3"2  times  in  the  duration  of  the  display.   In  the  ras- 
ter patterns  the  values  of  the  second  order  polynomial  mean 
square  fit  and  the  interpolated  raster  phase  for  one  bias 
condition  at  the  12  8  points  used  for  interpolation  are  shown 
as  the  characters  (+)  and  (*)  respectively. 

The  effect  of  the  quantization  interval  is  apparent 
in  all  raster  displays,  but  it  is  encouraging  to  note  that 
a  25  Hz  frequency  component  of  amplitude  equal  to  one-half 
the  quantization  interval  (Figure  2  0)  may  be  detected  even 
in  the  presence  of  a  transient  during  the  initial  half  second 
of  the  signal  sample.   This  detection  ability,  however,  is 
strongly  dependent  on  the  apparently  uniform  distribution  of 
transition  times  within  the  quantization  interval.   For  this 
deterministic  case,  as  the  raster  period  approaches  an  integer 
multiple  of  the  quantization  interval,  the  ability  of  the 
measurement  technique  to  isolate  a  signal  of  small  amplitude 
as  a  single  maxima  deteriorates. 

5 .   Laguerre  Polynomial  Coefficients 

In  the  attempt  to  characterize  the  apparent  transient 
behavior  of  some  of  the  raster  displays,  a  set  of  Laguerre 
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polynomial  coefficients  were  calculated  over  all  the  points 
z'(t-)  of  the  initial  portion  of  the  signal.   The  coefficients 
L   of  the  first  five  Laguerre  polynomials  were  obtained  by 
stepwise  numerical  integration  over  the  first  800  quantization 
intervals. 


Ln  =    E    z'(t.)  exp[-pt. ] (At.) 
u    t±<800     x  11 


L   =    E    z' (t.) [2pt.-l]  exp[-pt.] (At.) 
x    t.<800  11 


L~  =    E    z"  (t.)  [2(pt.)2-4pt.+l]  exp[-pt. ]  (At.) 
t-<800     ill  li 


L,  =    E    z' (t.) [4(pt.)3-[6pt.]2+6pt.~l]  exp[-pt. ] (At.) 
t.<800        J 


L.  =    E    z'  (t.)  [l-(pt.)4~(pt.)3+12pt.-8pt,+l]  exp[-pt.  ]  (At.) 

t.<800  i-  ii 

1  V-26 


p  =  80'      Ati  "  ti+1  ~  fci 

During  the  transient  period ,  the  aliasing  effect, 
which  occurs  if  the  nominal  raster  period  does  not  match 
the  actual  period,  is  especially  pronounced.   The  data  may 
very  well  be  undersampled  in  this  region,  thus  no  attempt  was 
made  to  compute  these  coefficients  independently  for  each 
bias  condition  and  the  phase  unwrapping  procedure  was  applied 
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to  all  data  points.   This  process,  however,  operated  on  the 
premise  that  the  phase  difference  from  transition  to  transi- 
tion is  not  more  than  one-half  the  nominal  raster  phase. 
Even  if  both  positive  and  negative  transitions  are  used  for 
phase  unwrapping,  this  premise  may  not  be  valid. 

If  the  initial  transition  times  predicted  by  the 
steady  state  period  and  knowledge  of  the  signal's  underlying 
data  modulation  can  be  determined,  a  more  accurate  represen- 
tation of  transient  phase  might  be  obtained.   Unfortunately, 
knowledge  of  this  data  modulation  was  not  available  for  the 
purpose  of  this  thesis. 

D.   COMPUTER  FILES  AND  PROGRAMS 

The  programs  used  to  produce  the  feature  measurements 
described  in  the  previous  section  were  designed  to  operate 
on  a  pre-production  model  of  the  Parameter  Encoder  which 
incorporated  in  its  hardware  a.  Hewlett-Packard  2100A  computer 
with  twin  disc  drives  and  moving  head  disc  operating  system. 
The  normal  graphics  display  and  list  device  was  a  Tektronix 
storage  scope. 

The  programs  discussed  in  this  section,  while  written 
in  HP  Fortran,  make  use  of  the  Fortran-callable  executive 
routines  available  under  the  disc  operating  system.   These 
routines  permit  program  overlays  and  disc  file  input/output. 
Similarly,  a  library  of  Fortran-callable  utility  routines 
developed  at  Electromagnetics  Systems  Laboratory  to  control 
graphics  in  put  and  output  were  used  to  create  and  erase  the 
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scope  displays,  temporarily  halt  computations  and  provide 
keyboard  control  of  the  measurement  process. 

The  signal  data  base  used  as  input  to  the  measurement 
process  was  stored  in  the  disc  file,  PARF.   This  file  con- 
tained a  signal  index  number,  the  automatically  measured 
clock  rate,  the  time  difference  between  successive  signal 
transitions  and   the  total  number  of  transition  times 
measured  (see  Table  1) .   Since  this  data  format  differed 
slightly  from  that  normally  accepted  by  the  standard  raster 
display  routine  RASP,  an  interface  program  FACE  was  written 
to  convert  the  data . 

Following  the  initial  display  of  the  raster  signal, 
process  control  could  be  transferred  to  a  program  FOURD 
whose  function  was  to  call  the  interpolation,  measurement, 
and  display  subroutines  and  to  store  the  measured  parameters 
in  the  disc  file  PAR.   The  overall  process  flow  of  FACE, 
RASF  and  FOURD  is  shown  in  Figures  23  and  24 . 

1 .   Measurement  Process 

The  first  subroutine  called  by  FOURD  is  RESDU  (see 
Figures  25  and  26  for  overall  process  flow  )  .   This  subrou- 
tine calculates  a  minimal  phase  difference  representation 
for  the  points  of  both  bias  conditions  outside  the  transient 
region.   The  points  thus  represented  are  used  by  the  subrou- 
tine CURV  to  calculate  the  coefficients  of  the  polynomial 
mean  square  fit.   They  are  also  used  to  calculate  the  signal's 
mean  bias  and  its  total  and  intrinsic  variance.   A  cubic 
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TABLE  I 
Disc  File  PARF 

Size:   1542  sectors  -  128  Integer  words  per  sector 
Sectors  0-5   Directory  -  768  words 


Word 


1-768 


Use 


Signal  number 


Comment 

Order  in  directory 
corresponds  to  order 
of  signal  files 


Sectors  6-1542   Signal  Files    2  Sectors  per  signal 


Word 
1-2 
3 

4  -  256 


Use 


Comme^ 


Nominal  period   Floating  point 

Number  of 
measurements 


Transition 
time 


Time  difference  between 
successive  transitions 
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FACE 


RASF 


NEW 
PERIOD 


READ  SIGNAL  NUMBER 


FIND  SIGNAL  FILE 


FORMAT 


FIND  FILE  TYPE 


FIND  TIME  BASE 
SCALES  START  AND  END  TIME 


± 


CALCULATE  RASTER 

POINTS 

DISPLAY   RASTER 


KEYBOARD 
CONTROL 


EXPAND 
SCALE 


SAVE 
DISPLAY 


MEASURE 


QUIT 


j. 


OVERLAY 

PROGRAM 

FOURD 


FIGURE    23.      Process   Flow,    Interface   and   Raster   Display 
Programs    FACE    and   RASF 
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CALL  RESDU 


INTERPOLATE 

FIT  POLYNOMIALS 

MEASURE  BIAS 

DISPLAY 

INTERPOLATION 

POLYNOMIAL 

FIT 


QUIT? 


.  .CALL-FQUR2  _ 

128  POINT 
FAST  FOURIER 
TRANSFORM 


I 


CALL  FORD 

CALCULATE  AND  DISPLAY 
FOURIER  MAGNITUDES 


il. 


QUIT' 


JL. 


STORE  PARAMETERS  IN 
FILE  PAR 


•;- 


OVERLAY   PROGRAM 
FACE 


FIGURE  24. 


Process  Flow  Measurement  and  Display 
Program  FOURD 
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SELECT 
ALTERNATE 
BIAS  POINTS 


SET  BIAS  CONDITION 
K=l 


INITIALIZE 
VARIANCE 


CALCULATE 


t,    ,?,  At3 


CALCULATE 

MINIMAL  PHASE 
DIFFERENCE 


CALL  CURV 

POLYNOMIAL  FIT 


DISPLAY 

POLYNOMIAL  FIT 


CALCULATE  INTERPOLATION 
FTX(K,N)  N=1-128 


DISPLAY 

INTERPOLATION 


SUBTRACT  POLYNOMIAL 
CONTRIBUTION  TO 
FTX(K,N) 


COMPUTE 

MEAN  BIAS 

VARIANCE 

LAGUERRE 

COEFFICIENTS 


START 


----© 

(see  Figure  26) 


1 


RETURN 


FIGURE    25. 


Flow  Diagram  Interpolation  and  Measurement 
Subroutine  RESDU 


7  8 


n  =  n+1 


n  =  n+1 


DEL  =  . 5*DELL 

PHH  =  PHL  +  DEL 


i  =  MAX 
J  =  2 


PHL  =  t.-nT 
X(l,l)  =  PHL 


PHH  =  t-nT 

DEL.  =  PHH  -  PHL 


UPDATE 
MEAN,  VARIANCE 


NO 


X(15J)  =  PHH 
X(2,J)  =  ti 

1  DELL  =  DEL 


if- 


PHL  =  PHH 


START  =  TsT 
PERIOD  =  T 
TRANSITIONS 
t 


MAX 


J  =  J+l 

— —TfT- — 


FIGURE  26.   Minimal  Phase  Difference  Algorithm 
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interpolation  scheme  CURF  then  provides  the  uniform  samples 
used  by  the  fast  Fourier  transform  subroutine.   RESDU  also 
calls  the  subroutine  LAGER,  which  performs  a  stepwise  numeri- 
cal integration  to  estimate  the  coefficients  of  the  first 
five  Laguerre  polynomials. 

The  two  other  subroutines  called  by  FOURD  are  F0UR2 , 
a  standard  fast  Fourier  transform,  and  FORD,  a  routine  v/hich 
calculates  the  average  square  magnitude  of  the  Fourier  com- 
ponents and  displays  the  results  on  the  Parameter  Encoder's 
storage  scope.   The  disc  file,  PAR,  which  contains  all  the 
above  feature  measurements,  is  organized  as  shown  in  Table  II . 
2 .   Statistics  and  Formatting  Routines 

The  analysis  of  the  features  measured  by  the  routines 
indicated  in  FOURD  makes  extensive  use  of  identification  and 
analysis  software  previously  developed  for  the  Parameter  En- 
coder.  These  procedures  use  a  yet  another  disc  file  structure, 
MASTF  (Table  III) ,  as  their  input  data  base.   Additional  refor- 
matting and  selection  of  the  parameters  was  necessary  to  allow 
the  use  of  this  file  since  it  accepts  a  maximum  of  50 
parameters . 

The  program  CSTAT  provided  class  statistics  for  the 
measurements  stored  in  MASTF.   In  addition  to  reformatting 
individual  signal  parameters,  and  storing  them  in  the  disc 
file  MASTF,  the  subroutine  TRNSG  performed  the  important 
function  of  interpolating  the  existent  Fourier  magnitudes 
to  estimate  the  square  magnitude  components  at  2  Hz  incre- 
ments.  This  was  necessary  since  the  resolution  of  each  signal 
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TABLE  II 
Disc  File  PAR 
Size:   200  Sectors   (128  Integer  Words  per  Sector) 
Sector  0    Directory 
Word  Use  Comment 


1-128 


Signal  Number 


Order  in  directory  corres- 
ponds to  order  of  files 


Sectors  1-199   Signal  Parameter  Files    2  Sectors  per  File 


Word 

1-2 

3-4 

5-8 

9-18 

19-28 

29-38 

39-40 
41-70 
71-72 
73-74 
75-78 
79-208 

209-256 


Comment 
Floating  Point  SIG(l) 
Floating  Point  SIG(2) 


Use 

Nominal  Period 

Signal  Resolution 

(Not  Used) 

Polynomial  Coefficients   Bias  Condition  I  SIG(5-9) 
Mean  Square  Fit 

Polynomial  Coefficients   Bias  Condition  2  SIG  (10-14) 
Mean  Square  Fit 

Polynomial  Coefficients   Floating  Point  SIG (15-13) 
Laguerre 

Mean  Bias  Estimate 

(Not  Used) 

Variance 

Intrinsic  Variance 

(Not  Used) 

Fourier  Magnitude 
Components 

(Not  Used) 


Floating  Point  SIG  (20) 

Floating  Point  SIG (35) 
Floating  Point  SIG (36) 

65  Terms  SIG(40-104) 
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TADLE  III 
Disc  File  MASTF 
Size:   1524  Sectors    128  Integer  Words  per  Sector 
Sector  0    File  Information 

Use 


Word 
1 
2 
3 
5 
6 

14 
15-64 

Signal  File 

Word 
1 

2-3 
4 

5-104 
105-154 
159-384 

Class  File 

Word 
1-2 
3 
4 

5-104 
105-204 
205-254 

255-304 
305-384 


Number  of  Sectors  in  MASTF 
Number  of  Words  in  Directory 
Number  of  Words  in  File 
Scale  Code  (-1  No  scaling) 
Variance  Breakpoint 
Number  of  Parameters 
Parameter  Numbers 


Use 
Signal  Number 
Identification  Number 
(Not  Used) 
Parameters 
Parameter  Weight 
(Not  Used) 


Use 

Identification 

Number  of  Entries 

Number  of  Parameters 

Class  Means 

Class  Variance 

Number  of  Parameters  Used 
for  Classification 

Parameter  Numbers 

(Not  Used) 


Comment 


50  Floating  Point  Words 
1  If  used,  0  If  Not  Usee 


Comment 


Floating  Point 
Floating  Point 
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was  variable.   Ground  truth  classification  of  each  signal 
had  been  previously  determined  through  the  use  of  parameters 
not  included  in  this  thesis,  and  was  stored  in  a  small  Disc 
File,  CSIGF. 

E.   ANALYSIS  OF  PARAMETERS 

From  the  basic  group  of  features  obtained  by  the  measure- 
ment process,  a  group  of  50  parameters  was  selected  for  further 
analysis.   The  analytic  techniques  used  were  those  readily 
available  to  the  user  of  Parameter  Encoder  identification 
software.   Several  iterations  of  measurement,  observation, 
and  feature  analysis  were  conducted  before  arriving  at  the 
parameters  indicated  in  Table  IV.  Earlier  measurements  had 
involved  the  use  of  fourth  order  polynomial  mean  square  fits 
to  the  entire  raster  display,  and  the  use  of  Fourier  magni- 
tude coefficients  of  raster  pattern   input  uncorrected  by 
the  polynomial  fit. 

These  earlier  techniques  appeared  to  provide  reasonable 
measurements  for  approximately  6  0  percent  of  the  signals 
considered,  but  failed  catastrophically  for  the  remaining 
signals.   The  major  cause  of  failure  appeared  to  be  the  ina- 
bility of  the  minimal  phase  difference  representation  to 
track  the  raster  pattern  through  initial  transients  and 
anomalous  transition  times.   The  final  measurement  process 
discussed  in  the  previous  section  was  applicable  to  about 
80  percent  of  the  signals  attempted  and  resulted  in  measure- 
ments for  the  final  training  set  of  74  signals  distributed 
among  nine  classes  as  shown  in  Table  V. 
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TABLE  IV 
Parameter  Identification 
Parameter  #         Identification 


901 

902 

903 

904 

905 

906 

907 

908 

909 

910 

911 

912 

913 

914 

915 

916 

917 

918 

919 

Nominal  Period 

Length  of  Signal  Resolution 

0   Order  Polynomial  Coefficient  Bias  Cond .  2 

1st  Order  Polynomial  Coefficient  Bias  Cond.  2 

2nd  Order  Polynomial  Coefficient  Bias  Cond.  2 

0   Order  Polynomial  Coefficient  Bias  Cond.  1 

1st  Order  Polynomial  Coefficient  Bias  Cond.  1 

2nd  Order  Polynomial  Coefficient  Bias  Cond.  1 

1st  Laguerre  Coefficient 

2nd  Laguerre  Coefficient 

3rd  Laguerre  Coefficient 

4th  Laguerre  Coefficient 

5th  Laguerre  Coefficient 

Mean  Bias 

Total  Variance 

Intrinsic  Variance 

Fraction  of  Total  Variance  Explained  by  Intrinsi 

Fourier  Power  Component   0  Hz 

Fourier  Power  Component   2  Hz 


950  Fourier  Power  Component  64  Hz 


84 


TABLE  V 
Training  Set  Signal  Distribution 
Class  1   (Not  Used  in  Distance  Ratio) 


1447 

1300 

875 

126 

171 

172 

775 

264 

Class 

2 

477 

758 

676 

Class 

3 

770 

183 

192 

1483 

1488 

693 

1354 

1367 

1498 

818 

306 

307 

319 

1525 

938 

Class 

4 

1441 

1444 

483 

1309 

771 

1327 

1337 

1338 

742 

315 

1399 

925 

Class 

5 

• 

127 

155 

766 

826 

Class 

6 

489 

235 

236 

787 

1543 

Class 

7 

1286 

1289 

672 

475 

674 

675 

120 

1304 

1305 

141 

Class 

8 

(Not  Used  in  Distance 

Ratio) 

888 

239 

1487 

698 

915 

446 

970 

Class 

9 

121         754         145        149         196 
189         201        202         203        204 
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The  ground  truth  identification  of  these  signals  was 
based  on  the  agreement  of  two  external  identification 
schemes.   One  of  these  was  the  semi-automated  process 
developed  by  Electromagnetic  Systems  Laboratory  which  uses 
additional  signal  parameters  not  associated  with  the  raster 
display.   The  other  process  uses  all-source  information  in 
arriving  at  signal  identifications. 

The  analysis  began  with  the  measurement  of  single  fea- 
ture separating  capability  for  each  of  the  50  selected 
parameters.   The  eleven  features  showing  the  greatest  capa- 
bility were  then  examined  for  redundant  features  by  calcula- 
tion of  class  and  global  correlation  matrices.   Finally,  the 
error  rate  of  a  minimum  distance  classifier  was  used  as  the 
criteria  for  a  feature  search  procedure  which  combined  single 
feature  ranking  and  search-without-replacement  techniques. 
Since  the  amount  of  training  data  was  limited,  even  this 
iterative  procedure  may  place  too  great  an  emphasis  on  the 
training  set  values. 

1.   Single  Feature  Separability 

The  single  feature  separability  measure  available 
in  the  Parameter  Encoder  software  is  a  modification  of  the 
distance  ratio  techniques  discussed  in  Appendix  C.2.   The 
program  FINTR  calculates  the  average  square  distance  between 
points  of  m  different  classes, 

K    K 

2     2    m-1  m     ,     i   j    (i)   (i)  2 

B     m(m  1)  .    j=.   K.  K  1 
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The  m  class  average  within-class-variance, 

m  Kk   Kk 

Dw  -  h     l     v    (v2   n   E   z   <xik)  -  x!k))2         V-28 
w   m  .  ,  K.  (K.  -1)  .  -  .  .    1       i 

k=l   k  k    i=l  3  =  1  J 


and  uses  the  ratio 

F  =  -S-j  V-29 

V 

as  a  measure  of  separability.   Since  the  array  storage 
available  in  this  program  limits  the  number  of  input  samples 
to  60,  the  signals  of  classes  1  and  8  were  not  used  in  the 
calculation  of  the  feature  distance  ratios.   The  distance 
ratios  obtained  using  all  signals  of  the  other  7  classes 
are  shown  in  Table  VI.  The  distance  ratios  of  the  frequency 
components  are  also  plotted  in  Figure  2  7. 

By  way  of  comparison  Figure  2  8  shows  the  distance  ratio 
as  a  function  of  frequency  for  data  obtained  in  earlier 
investigations  using  a  smaller  number  of  signals  and  classes. 
The  well  defined  maxima  in  the  latter  case  can  be  attributed 
to  the  fact  that  the  signals  and  classes  user  were  particu- 
larly noise  and  transient  free.   In  these  cases  one  was  able 
to  use  the  entire  raster  display  uncorrected  by  a  polynomial 
fit  to  estimate  the  Fourier  components. 

If  one  assumes  that  the  parameters  are  statistically  inde- 
pendent, then  the  ranking  of  features  according  to  distance 
ratio  provides  a  guide  for  feature  selection.   Eliminating 
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TABLE  VI 
Single  Distance  Ratios  for  Parameters   901  -  950 

Ratio         #       Ratio         #       Ratio 


901 

10.59 

921 

1.32 

941 

1.14 

902 

11.68 

922 

1.14 

942 

1.08 

903 

1.87 

923 

1.35 

943 

.94 

904 

1.07 

924 

1.25 

944 

.93 

905 

1.07 

925 

1.45 

945 

.93 

906 

1.17 

926 

1.81 

946 

.84 

907 

1.00 

927 

2.10 

947 

.95 

908 

1.03 

928 

1.51 

948 

.92 

909 

1.15 

929 

1.09 

949 

.96 

910 

1.03 

930 

1.44 

950 

.97 

911 

1.01 

931 

1.68 

912 

.99 

932 

1.06 

913 

1.08 

933 

1.27 

914 

.98 

9  34 

1.11 

915 

.99 

935 

1.28 

916 

1.20 

936 

1.44 

917 

1.14 

937 

1.35 

918 

1.68 

938 

1.26 

919 

1.22 

939 

1.02 

920 

1.19 

940 

1.01 
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FIGURE  27.   Feature  Separation  Quality  vs.  Frequency/ 
59  signals ,  7  classes 
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FIGURE  28.   Feature  Separation  Quality  vs.  Frequency, 
30  signals,  4  classes,  no  polynomial  fit 
corrections 
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those  features  whose  ratio  was  less  than  one  and  ranking 
the  best  eleven  of  the  remaining  resulted  in  the  feature  set 
shown  in  rank  order  in  Table  VII. 
2 .   Multiple  Feature  Analysis 

To  test  the  assumption  that  the  parameters  selected 
by  single  feature  ranking  were  independent  a  utility  routine 
PCORL  was  used  to  calculate  their  class  and  global  correla- 
tion coefficients.   Two  dimensional  scattergrams  for  selected 
pairs  of  highly  ranked  features  were  used  for  visual  confir- 
mation of  feature  correlation  and  clustering.   In  the  global 
case,  two  of  the  parameters  exhibited  correlation  coefficients 
in  excess  of  0.9.   These  parameters,  whose  scattergram  is 
shown  in  Figures  29  and  30,  are  the  nominal  raster  period  and 
the  signal  resolution.   Correlation  between  these  two 
parameters  might  be  expected  to  be  high  if  the  sample  signal 
length  corresponds  to  the  same  number  of  bauds  for  all  signals 
Both  the  scattergram  and  the  class  correlation  coefficients 
indicated  that  the  period  and  signal  duration  remain  highly 
correlated  on  a  class  basis.   Other  scattergrams  such  as  that 
shown  in  Figures  31  and  32,  for  the  next  most  effective  single 
feature  and  the  nominal  raster  period,  gave  little  evidence 
of  correlation.   The  results  of  these  calculations  and  obser- 
vations  for  other  feature  pairs  were  sufficiently  encouraging 
that  single-feature  distance  ratios  were  used  to  choose  the 
next  four  features  to  be  discarded  from  the  set  of  eleven. 
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TABLE  VII 
Ranked  Features 

J.  -  Ratio  Identification 

902  11.68  Signal  Resolution 
901           10.59            Period 

927  2.10  18  Hz  Fourier  Component 

903  1.83  0  Order  Poly  Mean  Square  Fit, 

Bias  1 

926  1.81  16  Hz  Fourier  Component 

931  1.68  26  Hz  Fourier  Component 

918  1.68  0  Hz  Fourier  Component 

92  8  1.51  2  0  Hz  Fourier  Component 

925  1.45  14  Hz  Fourier  Component 

936  1.44  36  Hz  Fourier  Component 

930  1.44  24  Hz  Fourier  Component 
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FIGURE  32.   Expanded  View  of  Figure  31 
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The  six  features  finally  chosen  were  selected  by 
repeated  application  of  a  minimum-normalized  distance 
classification  program  CNFUS .   Sixty-two  signals  from  all 
classes  except  class  four  were  classified  by  the  program 
according  to  their  distance  from  the  class  means.   This  dis- 
tance d.  of  a  particular  set  of  n  signal  features,  (x, . . .x  ) 
1  3  In 

from  the  class  mean   (y.   .  ..y    )  was  given  by 

In        3      J 

i  (i)v2 

2    !   n   (x   -  y    ) 

a.z  =  ±    i    — 3 — j-U—^ —  v-30 

1    n  j=l     (a?l})2 

.3 

The  signals  used  for  feature  selection,  their  distance 
to  the  nearest  four  class  means,  and  the  confusion  matrix 
resulting  from  the  choice  of  the  nearest  class  are  shown  in 
Tables  VIII  and  IX.   This  choice  of  six  features  leads  to 
approximately  56  percent  agreement  between  ground  truth  data 
and  the  first  choice  of  the  classifier  when  the  entire  train- 
ing set,  including  class  four,  is  used.   This,  of  course, 
represents  an  optimistic  estimate  of  the  probability  of 
correct  classification.   Further  testing  of  the  classifier 
over  an  independent  set  of  signal  data  is  advisable.   That 
the  addition  of  features  other  than  the  basic  clock  rate  does 
improve  classifier  performance  is  evident  in  the  improvement 
of  the  34  percent  classification  probability  obtained  using 
only  this  feature  as  the  basis  for  classification.   The 
resultant  single  feature  confusion  matrix  is  shown  in  Table 
X. 
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TABLE  VIII 

Classifier  Selections  of  4  Classes  Nearest 
to  Signals  of  the  Training  Set 


Parameters  Used 


902 

903 

927 

926 

931 

918 

Ground 

ignal 

Truth 

Classes  and 

Distance  to 

Mean 

Class 

427 

2 

2 

1.16 

1 

1.25 

4 

2.52 

5 

4.43 

758 

2 

2 

1.15 

1 

2.80 

676 

2 

1 

.70 

2 

1.12 

4 

1.96 

770 

3 

3 

.57 

8 

.73 

4 

.79 

6 

1.36 

183 

3 

4 

.51 

1 

.88 

8 

.94 

3 

1.07 

192 

3 

1 

.55 

3 

.84 

4 

.91 

8 

1.08 

1483 

3 

3 

.65 

4 

.66 

1 

.84 

8 

.85 

1488 

3 

8 

.74 

4 

.90 

3 

.99 

1 

1.23 

693 

3 

4 

1.07 

6 

1.23 

8 

1.44 

1 

1.52 

1354 

3 

4 

.72 

8 

.77 

3 

.78 

9 

1.08 

1367 

3 

3 

1.01 

4 

1.12 

1 

1.15 

8 

1.16 

1498 

3 

4 

1.24 

3 

1.25 

6 

1.46 

8 

1.49 

818 

3 

3 

.83 

4 

.95 

8 

1.23 

1 

1.28 

306 

3 

4 

.40 

3 

.62 

8 

.64 

6 

.82 

307 

3 

3 

.76 

4 

.86 

9 

.98 

8 

1.13 

319 

3 

3 

.96 

4 

1.14 

8 

1.32 

1 

1.78 

1525 

3 

8 

.57 

1 

.68 

4 

.81 

3 

1.08 

930 

3 

3 

.54 

8 

.56. 

4 

.84 

1 

1.05 

127 

5 

5 

.93 

1 

1.82 

4 

1.88 

8 

2.33 

155 

5 

5 

.83 

4 

1.40 

7 

1.51 

6 

1.71 

766 

5 

5 

.91 

6 

2.41 

4 

2.53 

8 

2.57 

826 

5 

8 

.42 

1 

.57 

4 

.78 

9 

.89 

489 

6 

1 

.84 

4 

.89 

6 

.92 

8 

1.17 

235 

6 

6 

.88 

4 

.99 

8 

1.61 

3 

1.65 

236 

6 

6 

.96 

4 

1.67 

1 

2.23 

5 

2.29 

98 


787 

6 

6 

1.21 

1 

1.33 

4 

1.56 

1543 

6 

8 

.70 

6 

.73 

4 

.85 

1 

.89 

1286 

7 

7 

.80 

4 

1.03 

1 

1.20 

6 

1.55 

1289 

7 

7 

.88 

4 

1.24 

1 

1.41 

5 

1.59 

672 

7 

7 

.88 

4 

.94 

1 

1.00 

9 

1.30 

475 

7 

7 

.61 

4 

1.38 

1 

1.59 

6 

1.76 

674 

7 

7 

1.25 

4 

1.34 

6 

1.51 

9 

1.62 

675 

7 

7 

.95 

1 

.96 

4 

1.57 

5 

1.05 

120 

7 

7 

1.37 

4 

1.59 

1 

1.91 

6 

2.28 

1384 

7 

.47 

4 

.79 

7 

.83 

6 

1.13 

1305 

7 

1 

.51 

7 

.65 

4 

1.2  3 

6 

1.92 

141 

7 

1 

.40 

7 

.99 

4 

1.14 

5 

1.53 

121 

9 

4 

.73 

9 

.89 

7 

1.05 

8 

1.14 

754 

9 

1 

.50 

4 

.53 

7 

.72 

5 

.94 

145 

9 

1 

.65 

7 

.82 

4 

.93 

9 

1.02 

149 

9 

3 

.75 

4 

.85 

8 

.87 

9 

1.02 

196 

9 

9 

.32 

3 

.59 

4 

.71 

1 

.74 

189 

9 

9 

.63 

4 

.87 

3 

1.08 

7 

1.09 

201 

9 

8 

.80 

1 

.80 

4 

.83 

9 

1.05 

202 

9 

4 

.56 

3 

.58 

8 

.58 

6 

.83 

203 

9 

9 

1.08 

4 

1.26 

1 

1.41 

3 

1.43 

204 

9 

1 

.93 

9 

1.00 

4 

1.25 

3 

1.32 

1447 

1 

1 

1.52 

4 

5.19 

2 

5.24 

1300 

1 

1 

1.52 

875 

1 

1 

.73 

4 

2.79 

2 

3.41 

8 

4.68 

126 

1 

1 

.40 

4 

3.75 

171 

1 

1 

.60 

4 

3.22 

2 

4.19 

172 

1 

1 

.79 

4 

3.29 

2 

3.41 

8 

5.64 

775 

1 

1 

.87 

4 

3.09 

2 

4.12 

264 

1 

1 

.62 

4 

2.56 

2 

2.67 

8 

4.37 

888 

8 

8 

.61 

1 

.83 

4 

1.30 

6 

1.79 

239 

8 

8 

.55 

3 

.69 

1 

.84 

4 

.89 

1487 

8 

8 

.76 

1 

1.16 

4 

1.24 

6 

1.58 

698 

8 

6 

.87 

4 

1.05 

8 

1.05 

1 

1.18 

915 

8 

8 

.90 

5 

.92 

1 

1.24 

4 

1.34 

446 

8 

1 

.66 

8 

1.04 

4 

1.12 

5 

1.29 

970 

8 

4 

1.00 

3 

1.02 

8 

1.32 

1 

1.34 
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TABLE  IX 

CONFUSION  MATRIX  FOR  CLASSIFIER  USING 
PARAMETERS  902,  903,  927,  926,  931,  918 


CLASSIFIER  CHOICE 
4  5  6  7 


10 


2 


3-1 
7 


GOOD   37  MISSED  25 

59%  correct  classification  over  training  set. 


3 
0 


#  missed 
Type  I 
Error 

0 
1 
8 
0 
1 
2 
3 
3 
7 
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TABLE    X 
CONFUSION    MATRIX    FOR   CLASSIFIER   USING   CLOCK    RATE    ONLY 

CLASSIFIER  CHOICE 

#  missed 
123456789  Type  I 

Error 
18--------0 


31-5  3---5110 

M       5        -  -  -  1  -  1  1  -  -        3 

w 
(d 

U61--3-1---4 

1 

4J7--12  12--4       10 

82-1--1-3-4 

9--  122-3119 

#of  503  11  34465 

excess 

Type  II 

Error 

#  GOOD   21      #  MISSED  41 

33.87%  correct  classification  over  training  set. 
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The  sequence  by  which  the  final  six  features  were 
obtained  involved  first  the  calculation  of  the  classifier 
confusion  matrix  as  each  of  the  four  lowest  ranked  features 
were  removed  from  the  set  of  eleven.   Their  removal  improved 
the  classifier's  performance  from  50  percent  to  58  percent. 
From  this  set  of  seven  features,  the  effect  of  removing  one 
parameter  and  using  the  remaining  six  yielded  the  following 
classifier  performance: 

Parameter  Removed       Probability  of  Correct  Classif icaticr 

58  % 
25  % 
53  % 

50  % 

51  % 
53  % 

59  % 
56  % 

The  above  suggests  that  the  best  single  feature  to  be  removed 

at  this  stage  is  9  01,  the  nominal  raster  period.   That  this 

remaining  group  of  six  features  is  not  necessarily  the 

optimal  choice,  can  be  shown  by  observing  that  a  choice  of 

five  features,  including  parameter  901  leads  to  a  61  percent 

probability  of  correct  identification  shown  in  Table  XI. 


None 

918 

931 

926 

903 

927 

901 

902 
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TABLE    XI 

CONFUSION    MATRIX    FOR   CLASSIFIER 
PARAMETERS    USED    901,    903 ,    926,    931,    918 


CLASSIFIER  CHOICE 

#  missed 
123456789  Type  I 

Error 
18--------0 


W5----3--1-1 
w 

"     6        1  -  -  -  -  4  -  -  -        1 


3 


92--3---146 

#  Excess    8009100600 

Type  II 

Error 


#  GOOD     38  #  MISSED     24 

61.29%  CORRECT  CLASSIFICATION  OVER  TRAINING  SET 
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VI. .  SUMMAPY  OF  CONCLUSIONS  AND  RECOMMENDATIONS 

It  has  not  been  the  intent  of  this  thesis  to  suggest 
feature  measurements  which  can  replace  those  currently  used 
in  the  Parameter  Encoder's  source  identification  scheme.   The 
fact  that  the  measurement  process  discussed  in  this  thesis 
is  currently  unable  to  provide  meaningful  measurements  for 
about  20%  of  the  signal  data  attempted  is  but  one  reason  for 
not  adopting  the  indicated  Fourier  components  and  polynomial 
coefficients  as  standard  features  to  be  measured. 

These  features  however,  constitute  a  hitherto  unexploited 
source  of  raster  pattern  information.   The  use  of  four  of 
these  features  in  conjunction  with  the  nominal  raster  period 
or  signal  clock  period  provides  approximately  100%  improve- 
ment in  classifier  capability  over  the  separation  performance 
of  the  nominal  raster  period  alone.   In  view  of  the  amount  of 
overlap  in  the  feature  distributions  for  several  classes, 
reflected  in  both  the  scattergrams  and  the  distances  from 
the  individual  signals  to  the  nearest  class  means,  it  is  not 
advisable  to  use  more  than  four  or  five  such  features.   Neither 
is  it  advisable  to  draw  too  firm  a  conclusion  as  to  the  opti- 
mal choice  of  these  new  features.   The  62  test  and  training 
signals  do  not  by  any  stretch  of  the  imagination  constitute 
a  completely  representative  statistical  base. 

It  is  rather  surprising  that  those  terms  intended  to 
characterize  the  mean  bias,  variance  measures  and  transient 
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phenomena  show  as  little  single- feature  separating  capability 
as  they  do,  since  it  was  precisely  these  phenomena  which  have 
been  used  by  system  operators  as  visual  aids  to  signal  iden- 
tification.  While  one  might  reasonably  associate  the  zero 
order  polynomial  mean  square  fit  coefficients  and  the  zero-Hz 
Fourier  component  with  transient  related  effects,  the  higher- 
order  Fourier  series  terms  have  not  previously  been  noted  as 
particularly  prominent  in  the  raster  signal  display.   In 
view  of  the  series  of  nonlinearities  present  in  obtaining  the 
interpolated  minimal  phase  raster  points,  it  is  not  impossible 
that  these  higher  order  terms  owe  their  separating  capability 
to  the  data  modulation  of  the  signal  source.   In  this  respect 
the  "sideband"  structure  of  the  single  feature  separation 
data  of  Figure  27  is  most  intriguing. 

Before  completely  discounting  the  usefulness  of  bias 
measures  and  transient  phenomena,  a  more  careful  represen- 
tation and  analysis  of  the  data  in  the  transient  region  should 
be  attempted.   The  representation  could  perhaps  involve  the 
calculation  of  the  minimal  phase  by  use  of  a  moving-average 
technique  which  works  backward  from  data  in  the  steady-state 
region  to  the  starting  time  of  the  signal.  At   the  same  time 
use  of  a  different  decay  constant  and  a  smoother  integration 
to  obtain  the  Laguerre  polynomial  coefficients  might  also 
prove  of  value. 

If  additional  accurate  ground  truth  data  is  available, 
the  modeling  of  the  probability  distributions  by  something 
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other  than  a  Gaussian  form  could  be  useful.   The  average 
squared  magnitude  of  the  Fourier  components  might  be  more 
accurately  represented  by  a  Rayleigh  distribution.   Although 
some  of  the  measured  parameters  for  one  or  more  classes  may 
not  have  a  unimodal  distribution,  the  number  of  accurately 
labeled  samples  currently  available  provides  little  encour- 
agement for  the  use  of  Parzen  window  probability  estimators. 

The  foregoing  comments  suggest  that  four  or  five  parameters 
characterizing  the  raster  scan  data  may  provide  signal  iden- 
tifying information  useful  to  about  the  same  degree  of  accura- 
cy as  visual  interpretation  of  the  raster  scan  display.   Since 
the  optimal  parameters  measured  in  this  investigation  do  not 
appear  to  correspond  to  the  visual  cues  used  by  system  opera- 
tors, one  has  reason  to  believe  that  other  parameters 
derived  in  subsequent  analysis  might  improve  upon  this 
degree  of  accuracy. 

The  measurement  technique  investigated  exists  outside  the 
mainstream  of  the  Parameter  Encoder's  indentif ication  process. 
As  such  it  constitutes  a  relatively  slow  procedure  (one  minute 
per  signal) ,  requiring  the  use  of  disc  data  files  which  con- 
tain much  of  the  same  information  available  directly  from 
the  master  file  structure.   Even  should  the  programs  used 
be  modified  to  work  solely  with  the  master  file,  the  data 
obtained  might  be  of  greatest  use  in  a  sequential  identifica- 
tion scheme.   It  could  be  used  much  as  the  present  raster 
scan  display  is;  to  resolve  ambiguities  in  the  choice  between 
a  relatively  small  number  of  classes. 
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APPENDIX  A 
METHODS  OF  PROBABILITY  ESTIMATION 

1 .   Parametric  Techniques 

Perhaps  the  most  common  technique  which  has  been  used 
to  infer  the  a  priori  probability  structure  presupposes  that 
the  distribution  of  features  in  a  given  class  is  multivariate 
normal.   As  a  first  step  the  mean  and  variance  for  each  feature 
of  a  class  are  estimated.   For  the  class  C.  having  samples 
(x,  ...  yL.    )  ,  these  estimates  may  be: 

„<«  =  1   ^  s(i)  A-! 

-       Ki  j=i  -^ 


K. 

(i)2  _   1    J1  ,        (i).2  _  0 

ak    =  K7^T  *      (xjk  "  ^k   )  A"2 

i    3=1    J 


(i)    ,   (i)2  c      i      2(i)    mU 

where  y,    and  a,     are  components  of  y   and  o  .   The 

K  K  —        — 

variance  vector  represents  a  simplification  from  the  more 
general  multivariate  form: 


p(x|C.)  =  — -yl exp[-y(x  -  y1)TsT1(x  -  ufl))l     A-3 

2-rry/  IS.  I  -11 

The  components  of  the  vector  a     are  the  diagonal  elements 
of  the  d  x  d  covariance  matrix,  which  may  be  estimated  by: 

K. 

(i)   ..(i)w„i   ..UKT         A_4 


S.  =  _i_  z   (xUJ  -  yU;)(x1  -  y!1') 
~l    K.-l  •  ,   ~--j      —1     —1  —j 

l    j=l    J       j      j     j 
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Under  the  assumption  of  independent  features,  the 
distribution  is  characterized  by  the  product  of  d  Gaussian 
distributions . 

d 

p(2£|c.)  =   n  T expf-l  -^  (xk-y£)2]       A-5 

K~±  ,9  J  (i)        ak 

If  required,  the  technique  may  be  extended  to  include  the 
estimation  of  a  covariance  matrix  for  each  class  by  use  of 
Equations  A-3  and  A-4 .   This  procedure  has  enjoyed  great 
popularity,  particularly  in  the  modeling  of  features  with 
unimodal  class  conditional  probability  distributions.   The 
ease  with  which  the  probability  measure  of  a  Gaussian  distri' 
bution  may  be  interpreted  as  a  weighted  distance  measure 
has  been  an  important  factor  in  its  retention. 
2 .   Non-Parametric  Techniques 

Another  technique  technique  for  obtaining  the  a- 
priori  probability  density  function  involves  its  direct 
estimation  through  the  use  of  window  functions  or  potential 
functions.   In  this  technique,  the  estimate  at  some  point  x 
in  the  d-dimensional  feature  space  is  obtained  by  the  super- 
position of  terms  which  are  generalized  distance  functions 
of  that  point  and  the  training  set  points  y .    of  the  class 


C.  . 

1 


K. 

1    1       ^ 

p(x|ci)  =  y~      E  Y(X/^-;)  a"6 

i  j=l       J 
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The  technique  arises  from  a  viewpoint  proposed  by 
Parzen.   In  one-dimensional  cases,  one  can  estimate  the 
probability  density  at  some  point  x  by  counting  the  number 
of  samples  N  (x)  in  the  interval  [x+h,x+h]  and  dividing 
by  the  product  of  the  total  number  of  samples  M  and  the 
interval  width  2h 


pM(x)  -  -W  A"7 


If  one  defines  a  "window"  function 


k>  |y|  i1 

K(y)  =    z  A-8 

0   lyl  >  1 


one  can  write 


,   M    x-y . 


This  corresponds  to  the  potential  function 


Y<x,y)  =  h  K  (2L-^)  A"10 


Additionally,  the  argument  presented  may  be  extended  to 
multiple  dimensions  by  considering 
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1   |u. I  <  T  for  j  =  1  ,...,m 
K(u)  =     '   3   ~  l 

0  >    otherwise 


i        2£  -  Z 

(x,y)  =  ~  K  (-^ )  A(ll-12) 

m      m 


where  V   is  a  multidimensional  volume  of  the  hypercube  of 

m  J  c 

side  h  . 
m 

The  window  function  K  need  not  be  the  exact  form  of 

Equation  A-ll.   It  can  be  shown  that  if  the  estimate  PM(x) 

is  a  function  both  of  the  total  number  of  samples  M  and 

V  =  V  (M) ,  the  conditions 
m    m 


0  <  K(y)  <  °° 


/K(y)  dy  =  1 


lim  V  =  0 


lim  MV   =  °° 


A(13-16) 


ensure  that  as  the  number  of  total  samples  M-*-«> 


lim  j3M(x)  =  p(x)  A-17 

fy\-yco 


lim  (pM(x)  -  p(x))2  =  0  A-18 

M-)-oo 


110 


Thus,  we  can  consider  PM(x)  to  be  a  "blurred"  or  "noisy" 
version  of  p (x)  as  seen  through  an  averaging  window. 

Intuitively,  one  can  describe  the  following  useful 
properties  of  the  potential  function:   y(x,y) 

1.   y(x_/Y)  should  be  maximum  for  x  =  v_. 
•  2.   y(x_/Y.)  should  be  approximately  zero  for  x  distant 
from  y. 

3.  Y(x,y)  should  be  continuous  and  decrease  monotonically 
with  the  distance  between,  y  and  x. 

4.  If  y(2£i  »y)'  =  Y^o'^J  w^ere  y  is  a  sample  point,  the 
patterns  represented  by  x,  and  x~  should  have  about 
the  same  similarity  to  y. 

Such  properties  can  be  obtained  with  a  potential  function 
of  the  form: 

1  1    a     (x.  -  yi)2 

Y(x,y_)  =  — ^72 exp[-  ^      z 2 3     A~19 

2iT/a....a,  i=l    a  • 

Id  i 

where  a.  are  arbitrarily  chosen  values.   Although  the  form 
here  is  Gaussian,  and  the  function  of  a  distance  measure, 
it  is  also  possible  to  construct  a  potential  function  of 
the  form: 


N     2 

Y(x,y)  =  I  X/  <j).  (x)  <J>,  (y)  a-20 

i=l  x   x     x 


where  {$.}  is  a  complete  set  of  multivariate  orthonormal 
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functions  and  A.  are  constants.   Such  techniques  have  been 
described  more  completely  in  the  work  of  Aizermann  and 
Bravermann  and  others  [1,  7,  19]. 

Given  sufficient  samples,  the  Parzen  window  or 
potential  function  approach  essentially  assures  satisfactory 
convergence  to  an  arbitrarily  complex  distribution.   Unfor- 
tunately, this  sufficient  number  of  samples  may  be  far  greater 
than  the  number  required  if  the  form  of  the  distributions 
was  known.   Since  every  sample  point  is  used  in  the  construc- 
tion of  the  density  function,  the  above  approach  affords 
little  economy  in  the  way  of  data  reduction  and  leads  to 
a  demand  for  computation  time  and  storage  space  exponentially 
increasing  with  the  number  of  features. 
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APPENDIX  B 
CLASSIFICATION  ALGORITHM 
1 .   Distance  Based  Classification  Schemes 

An  assumption  commonly  used  to  implement  a  decision 
rule  is  that  of  multivariate  normal  class  dependent  density 
functions  of  the  form: 


p(x|C.)  =  p(C.)— =rJ exp[-  y[x-y  (l)  ]  V"1  [x-y  (l)  ] 

2  TT  7   I  S  .  I  ^  X   - 


=  p(Ci)  ~-  exp  [-  |[x-y(l)]TSi1[x-y(l)]    B-l 
i 


where 


x  is  the  d-component  feature  vector 
y   is  the  d-component  class  mean  vector 
S.  is  a  d  x  d  class  covariance  matrix 
|S. I  is  the  determinant  of  S.. 

I  i  ■  — J. 

For  such  a  distribution  one  may  obtain  the  locus  of 

points  of  constant  density  as  a  hyperellipsoid  for  which 

(i)  T  -1     (i) 
the  quadratic  surface  (x-y    )  S.  (x-y_   )  is  constant. 

The  quantity 


2    ,    (i).T  -1,    (i)  .  n  9 

r.   =  (x-y   )  S^    (x-y_   )  B-2 
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is  sometimes  called  the  Mahalanobis  distance  from  x  to  y 

Since  the  probability  densities  are  always  non- 
negative,  the  maximum  likelihood  criteria 


P(Ci)p(Ci|x)  >  p(Ci)p(C.  |x)    for  all  j  B-3 

implies 

log  p(Ci)p(Ci|x)  >  log  p  (C^p  (C^  |.x)  B-4 

Substituting  the  multivariate  density  from  Equation 
1,  the  above  expression  yields 


log  p(Ci)  +  log  V±  -  ^[x-/]TS  1[x-y1] 


*  !£J  (log  P(C.)  -  log  V.  -  |[x-l(i)lT£j1trii(i)l> 

B-5 


The  above  equation  defines  the  decision  boundaries 
of  each  class.   If  a  common  covariance  matrix  S.  =  S .  =  S 
is  assumed  for  each  class  in  the  above  equation,  the  log  V. 
terms  may  be  ignored  and  the  form  of  the  decision  boundary 
becomes 


1,    (iKTc-l,    (i). 
^-(x-y    )  S       (x-y    ) 

<   min{  [x-yu;]  S   [x-yU;]  +  log  — }       B-6 

j^i   ~-         --  PtCj) 
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Similarly,  if  equiprobable ,  a-priori  class  proba- 
bilities are  assumed,  p(C.)  =  p(C-)  then  the  decision  rule 
becomes 


11)   [x-y1]TS~1[x-y(l)]    min{[x-y(j)]S~1[x-v(j)]}   B-7 


j*i 


This  is  equivalent  to  simply  assigning  the  feature  vector 
x  to  that  class  whose  pattern,  y    ,  is  at  the  minimum 
Mahalanobis  distance  from  x.   The  concept  of  Mahalanobis 
distance  is  particularly  applicable  to  well-separated 
unimodal  feature  distributions,  but  presents  a  considerable 
problem  in  those  cases  where  the  distributions  exhibit 
more  than  one  local  maximum. 

2 .   K  -  Nearest  Neighbor  Algorithms 

The  concept  of  distance  may  be  applied  in  a  slightly 
different  form  to  probability  functions  not  well  characterized 
by  the  multivariate  normal  form.   In  this  application,  dis- 
tances in  feature  space  are  quite  frequently  normalized  by 
some  form  of  global  variance  or  maximal  spread.   The  nor- 
malized Euclidean  distance  from  the  sample  point  x  to  each 

-  +_.  .  .        ,  (1)     1       (m)     (m)  .  . 

point  of  the  training  set  (x,   .  . .  xT.  ...,x,   .  .  .x„   )  is 

1  m 

calculated  and  ranked  in  order  of  increasing  distance.   In 

the  simplest  algorithm,  class  membership  of  the  sample  point 

is  determined  by  the  class  membership  of  the  greatest  number 

of  the  k-nearest  neighbors  in  the  training  set.   For  example, 

suppose  that  of  the  ten  points  of  the  training  set  nearest 
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to  the  test  sample,  three  were  members  in  class  1,  four  of 
class  2  and  one  each  were  members  of  classes  4,  5,  6,  the 
simplest  algorithm  would  choose  class  2  as  the  probable 
classification.   Alternative  algorithms  of  this  type  weight 
the  contribution  of  the  k-nearest  training  points  according 
to  their  distance  ranking.   The  nearest  point  might  be 
assigned  a  weighted  class  vote  of  10,  the  next  nearest,  a 
weight  of  9  and  so  forth.   Class  membership  of  the  sample 
point  would  then  be  determined  by  summing  the  weighted  class 
votes  of  the  ten  nearest  neighbors  and  choosing  that  class 
which  received  the  maximum  number  of  votes. 

Such  procedures  are  reminiscent  of  the  Parzen  window 
estimators  described  in  Appendix  A. 2,  and  in  fact  can  be 
shown  to  be  equivalent  if,  instead  of  rank  order,  distance 
weighted  contributions  to  the  class  vote  are  allowed.   Numer- 
ous techniques  [3,  4,  10]  exist  to  increase  the  computational 
efficiency  of  these  K-nearest  neighbor  algorithms  by  reducing 
the  number  of  Euclidean  distance  terms  which  must  be 
calculated  for  each  sample  point.   It  must  be  mentioned 
that  these  techniques  are  of  greatest  use  when  the  individual 
classes  are  equiprobable  or  at  the  very  least,  when  every 
class  of  the  training  set  has  more  than  k  elements. 
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APPENDIX  C 
FEATURE  SELECTION  AND  RANKING 

1 .   Linear  Combinations  of  Features 

The  objectives  of  low  dimensionality  and  retention 
of  sufficient  information  may  often  be  obtained  through  the 
use  of  custom  orthnormal  transformations  of  the  measurement 
space  [13,  15,  19,  20].   In  such  an  interpretation  the  sample 
data  is  presumed  to  be  adequately  represented  by  the  linear 
combination  of  some  finite  set  of  functions  _<£>,  .  .  ._£ 

n 
z  =   £  a.  d>.  C-l 

-   i=l  X  ^ 

where 

d).  =  (d).n  ,  .  #  #  i  it.  ) 

For  a  given  set  of  functions  the  weights  are  adjusted  so 
that  the  average  mean  square  error 


n         n        9 

I      (zn  -   E   a.<f> ik)Z  C-2 

k=l   n   i=l  X   1K 


is  minimized. 

The  values  a.  are  thus  determined  by  taking  partial 
derivatives  with  respect  to  a.  and  setting  the  result  equal 
to  zero,  obtaining  n  equations  in  n  unknowns  of  the  form 
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n      m  m 

1     ai  l      *ik  ^ik  =  Z  zk  *ik  C"3 

i=l  1  k-1      JK       k=l  K  1K 


Now  if  the  functions  _£....._£  are  orthonormal;  that  is  if 


m  0     i  i-   j 

Z  *ik  ^ik  =  C"4 

k=l  1K  3K  1    i  =  j 


Equation  C-3  reduces  to 


m 


a .  =   Z   z,  A. , 
1   k-1   k   ik 


C-5 


Then  values  of  a.  may  then  be  used  as  the  pattern  space  x 
representation  of  the  vector  z  under  the  linear  transformation 


x  =  T  z  C-6 

T  =  [jj>.  , ^2  '  •  •  •  'j£.n^ 

Expansion  of  a  function  by  a  Fourier  series  represents  one 

form  of  orthonormal  transformation,  but  one  which  has 

distinct  advantages  in  computational  efficiency. 

It  can  be  shown  that  Karhunen-Louve  expansion  of 

N  m-dimensional  samples  z.  =  (z.,,...,z.  )  leads  to  a  m  x  m 

r     l     ll      lm 

eigenvalue  problem. 


«  *k  -  °k  *k  c'7 
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where 


1  N 
krs  =  N  .^  zir  Zis  C"8 


£  "  [krsl 


The  above  solution  does  in  fact  minimize  mean  square  error 
averaged  over  the  entire  sample  set.   By  choosing  eigen- 
vectors corresponding  to  the  n  largest  eigenvalues,  one  can 
achieve  dimensionality  reduction  from  m  to  n.   This  approach, 
also  known  as  principal  factor  analysis,  requires  the  compu- 
tation of  the  eigenvectors  and  eigenvalues  of  a  large  (m  x  m) 
matrix.   Note  that  if  the  number  of  signals  N  is  less  than 
m,  there  will  be  at  most  (N-l)  non-zero  eigenvalues,  and  the 
minimum  mean  square  error  is  zero. 

A  simpler  procedure  is  to  choose  a  small  number  of 
linearly  independent  patterns  typical  of  the  different  classes 
or  types  of  samples;  then  apply  the  Gram-Schmidt  orthonormali- 
zation  process  to  them  directly.   In  this  process  one  begins 
with  the  "typical"  samples,  perhaps  derived  from  class  means. 


(1)    ,  (1)      (IK 

y    =  (y,   ...  y '  ) 
—        1        m 


(2)    ,  (2)       (2), 

y   =  (y,   ...  y'  ) 

—        1        m 


,  (n)  _  ,  (n)       (n) .  Q 

y    =  (v1         . . .  ym   )  C-9 
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and  constructs  the  orthonormal  functions  $.   in  the  following 


manner: 


where 


_1 

£1  =  %'!]_)   ^ 


(2) 
C12  =  ^l'ii  '   ) 


i2  -  h      "  ci2  ii 


c,  =  (A.,u(k)) 

j  k    x j  — 


4  =  u   "  cik^i  "  c2ki2---  "  ck-i  k  4-i 

_i 
lk  =  (VV  2  ^k  c"10 


i    m 

(g-  »9iJ  =  —  £  g.  .  ,g.  .  C-ll 

— 3  — k    m  •  _■,  ji  ki 

Clearly,  all  of  the  samples  used  to  generate  the  orthonormal 
vectors  may  be  represented  exactly  by  an  n-term  expansion; 
thus  the  representation  error  of  the  'typical'  sample  is 
minimized.   If  these  samples  are  truly  typical,  one  can 
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expect  other  samples  to  be  represented  closely  by  the 
orthonormal  basis  functions. 

The  dimensionality  of  the  resultant  pattern  space 
is  the  same  as  the  number  of  linearily  independent  typical 
samples . 

2 .   Criteria  for  Feature  Selection  and  Ranking 

Rather  than  achieving  a  pattern  space  through  use 
of  a  transformation  which  minimizes  mean  square  error  over 
the  training  set,  it  is  often  desirable  to  optimize  some 
measure  of  separation  quality.   Meisel  [14]  catalogs  more 
than  ten  different  measures  of  separation  quality  which  are 
directly  calculable  from  general  a  priori  conditional  density 
functions.   For  unimodal  distributions,  however,  it  is 
convenient  to  describe  such  a  measure  in  terms  of  average 
distances  in  the  feature  space. 

Two  distance  measures  need  to  be  considered.   One, 
the  inter-class  distance,  might  be  characterized  by  either 
the  distance  beteen  the  means  of  the  classes  or  the  average 
squared  distances  between  points  of  different  classes. 
Another,  the  intra-class  distance,  is  the  measure  of  the 
internal  scatter  of  samples  and  may  be  characterized  by  an 
average  squared  distance  between  the  points  of  each  class 
and  summed  over  the  classes.   In  the  two-class  problem,  a 
common  definition  of  the  interclass  distance,  S,  and  intra- 
class  distance,  R. ,  is 

1 
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K    K 

1  2       11)  (2) 

S  =  =~V-   Z    Z   d'[x.U',  x)Z)]  C-12 

Kl  K2  i=l  j=l     _1     _:I 

K    K 

2      -1   -,1  „2.  (i)    (i)  ,  _  ,  „ 

Ri  =  K.  K.  .    *    *  d  '2j   •  2k  ']  C-13 

i   1-1   3=1  k=3      J 

2 

where  d  (x,y)  is  some  distance  measure  of  the  vector 

variables   x,y  .  . 

Conventional  distance  measures  might  use  one  of 
the  following  forms: 


2         n  2 

Euclidian   d,  (x,y_)  =   Z  (x.  -  y.) 

1         i=l   1    1 


n 
city  block  d2 (x,y)  =  Z  |x.  -  y. | 

i=l 

2  V1 

Mahalanobis  d_  (x,y)  =  (x  -  y)  ~   (x  -  y) 


localized   d  *(x,y)  =  1  -  exp [~   d  2 (x,y) ]         C(14-17) 
distance 


In  generalizing  this  distance  criteria  to  a  m-class  problem, 
it  is  convenient  to  resort  to  the  concept  of  scatter  matrices 
where 

K. 

1         i         (i)  T 
Si  =   Z   (xk  -  y  ) (xk  -  y/  V  C-18 

k=l 


m 

S   =   Z  S.  C-19 

~w    .  . ~i 
1=1 
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represent  the  individual  class  scatter  matrices  and  the 
collective  within-class  scatter  matrices  respectively.   If 
one  defines  a  total  scatter  matrix  S   in  terms  of  a  global 
mean  M 


.    m       , . . 
M  =  ——"   E   K.  yU;  C-20 

i    X  1  =  1 


S   =   E  (x  -  M) (x  -  M)T  C-21 

all  samples 


we  find  that 


m  (i)        (i)  T    m  (i)  (i)  T 

S   =   E    E(x  -  yUM  (x  "  P   M   +  I         Z(Vl  '  -  M)  (U.U;  -  M)1 

i=l  xeC.  i=l  xeC 


!„  +       Sb  c"22 


provides  a  natural  definition  for  both  inter-set  and  intra- 
set  scatter  matrices: 


m 
S   =  I      S.  C-23 

~W    i=l   x 


m       (i)        i     T 
SD  =   E   K.  [\iK    '    -  M)  (y   -  M)1  C-24 

~B    i==1   l  -      -  _ 


These  three  matrices  not  only  provide  the  basis  for  multiple 
discriminant  analysis  [7,  14,  18],  but  their  trace  or 
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determinant  can  also  serve  as  a  measure  of  separability  for 
use  in  an  algorithm  involving  either  linear  combinations 
of  features  or  discrete  feature  selection.   For  example, 
as  criteria  for  separability,  one  might  choose  to  maximize 
the  quality  function 


Q.  =  tr  S   :  SD  C-25 

1       ~w    ~B 


or  to  minimize  one  of  the  quality  functions 

Q2  =  tr  ?T~1  5w  C"26 

|S  I 
Q,  =  — ~  C-27 

T 

All  of  the  above  criteria  are  particularly  interesting  in 
that  they  can  be  evaluated  as  a  result  of  multiple  discrim- 
inant analysis  using  linear  combinations  of  features. 

For  a  c-class  problem  we  can  accomplish  the  pro- 
jection of  the  d-dimensional  feature  space  onto  some  alterna- 
tive space  of  not  more  than  c-1  dimensions  through  the  use 
of  a  non-unique  d  x  (c-1)  matrix  W. 


y  =  WT  x  C-28 


It  is  straightforward  problem  to  show  that  the  resultant 
scatter  matrices  S  ,  S   in  this  projection  space  are 

~W   ~  13 
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s  =  wT  s  w 


S0  =  WT  SD  W  C-29 


If  one  then  chooses  W  to  maximize  the  ratio  of  between-to 
within-class  scatter  matrix  determinants, 


|sB|  |wT  S  W| 
Q4  =  -^-  =   ~   —  -  C-30 

Is  |w  s  wl 

1 ~w '  !.   ^w  ~ ' 


there  is  a  solution  of  the  form 


SR  w.  =  A.   S   w. .  C-31 

^.B  —1     l   ~w  —l 


If  S   is  non-singular  one  can  directly  solve  an  eigenvalue 
~w 

problem  related  in  form  to  the  Karhunen-Louve  problem 
Equation  C-7. 


(StT  1  S_)  w.  =  A.  w.  C-32 

v^w    ~B   —l     l  —l 


or  one  may  then  solve  the  characteristic  polynomial 


S„  -  A.  S  I  =0  C-33 


for  the  largest  eigenvalues  and  find  the  corresponding 
eigenvectors  from 
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<fB  "  Xi  !w»  2i  -  0  C-34 


If  m-1  of  the  vectors  (y    -  M)  are  linearly 
independent,  there  will  be  not  more  than  m-1  non-zero  eigen- 
values for  this  problem.   If  the  within-class  scatter 
matrix  is  isotropic,  which  is  the  case  if  all  features  are 
independent  and  suitably  normalized,  the  resultant  eigen- 
vectors span  the  space  defined  by  the  vector  set  (y    -  M) 
and  may  be  established  by  Gram-Schmidt  orthogonalization . 
The  eigenvalues  of  this  problem  serve  to  characterize  the 
trace  and  determinant  criteria.   For  example,  the  use  of 
all  features  leads  to  the  criteria  values 

,      m-1 

Q.  =  tr  S     S   =   Z    A. 
1        «w    ~B    .=1    1 

m-1     n 

1=1         1 

Is  I     m~1  i 

|St|      1=1  1 

As  was  previously  mentioned,  the  criteria  described 
here  are  best  suited  to  unimodal  distributions  and  increase 
the  theoretically  achievable  error  rate.  In  those  cases 
where  it  can  be  established  that  a  unimodal  distribution  is 
inappropriate,  the  use  of  separability  measures  more  closely 
related  to  the  estimated  conditional  probability  density  is 
indicated  [16] .  For  example  one  might  use  a  quality  measure: 
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-,      N    9  m  ~ 

-  i   I   [fT(F(x.))  -   Z   [p(Ci)  p(F(x,) |C±) 3  C-38 

i=l        J     i=l  J 


where  N  is  the  total  number  of  training  set  samples  and 
F(x)  represents  some  selection  of  features  from  the  initial 
feature  space  and 


m 
P(x)  =   E   p(C.)  p(x|c.)  C-39 

i=l  ± 


Qr  is  in  fact  a  measure  of  overlap  and  has  an  optimal  value 
of  zero. 

3.   Feature  Search  Algorithms 

A  cursory  inspection  of  the  problem  of  finding  those 
n  of  m  features  which  optimize  one  of  the  criteria  given  in 
the  previous  section  shows  that 


N<_  =  —r^ rr  C-4  0 

t   n! (m-n) ! 


evaluations  must  be  performed  if  one  is  to  conduct  an 
exhaustive  search  of  all  possible  feature  combinations.   The 
exhaustive  search  for  the  best  ten  or  twenty  features  for 
example,  requires  consideration  of  more  than  184,000  possibil- 
ities.  Typically,  the  practical  feature  search  employs 
suboptimal  procedures  that  may  be  justified  if  simplifying 
assumptions  regarding  the  nature  of  the  feature  distributions 
are  allowed. 
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The  simplest  suboptimal  search  procedure  is  to  assume 
the  independence  of  features  and  evaluate  the  desired 
criteria  function  for  all  single  features.   These  m  features 
are  then  ranked  in  order  of  decreasing  optimality  and  the 
first  n  are  chosen.   The  number  of  criteria  evaluations  is 
the  same  as  the  dimension  of  initial  feature  space. 

A  second  search  method  used  by  Mucciardi  and  Gose 

[17]  involves  a  technique  known  as  search  without  replacement 

This  method  requires  N_  criteria  evaluations  where 

hi 


Ne  =  nOnj-J}!  c_4i 


Fu  [9]  has  proposed  a  technique  of  sequential  feature 
selection  which  can  be  modified  to  include  the  cost  of 
feature  measurement  and  Chang  [3]  has  developed  alternative 
dynamic  programming  approaches  to  feature  selection,  one  of 
which  requires  only  a  slightly  greater  number  of  evaluations 
N„  than  the  search  without  replacement  algorithm. 


NE  =  n(m  -  ^^2^)  +  n  -  2  (m  -  1)  C-42) 


The  single  feature  ranking  technique  has  the 
disadvantage  of  ignoring  the  effects  of  feature  correlations. 
The  search  without  replacement  technique  provides  a  method 
for  treating  these  correlations,  but  it  also  presumes  that 
the  features  obtained  through  n  stages  of  conditional  single 
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feature  evaluation  procedure  are  the  n  best  features,  which 
need  not  be  the  case.   In  general,  as  one  progresses  to 
successively  more  complicated  search  procedures,  more  of 
the  effects  of  feature  dependence  may  be  taken  into  account 
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PARAMETER   HISTOGRAMS 
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COMPUTER  PROGRAMS 

THE  PROGRAMS  DISCUSSED  IN  THIS  SECTION  ALL  USE 
LIBRARY  AND  EXECUTIVE  ROUTINES  WHICH  ARE  PART  OF  THE  PARA- 
METER ENCODER  SOFTWARE  PACKAGE.   AS  AN  AID  TO  INTERPRETATION 
A  BRIEF  DESCRIPTION  OF  SOME  OF  THESE  ROUTINES  FOLLOWS 

ROUTINE        FUNCTION 

ALPHA  CONTROLS  TYPE  OF  OUTPUT  TO  CRT 

CHOUT(I)       OUTPUTS  SINGLE  CHARACTER  TO  CRT  BUFFER 

CURSR(  IC,IX,IY) 

DISPLAY  CURSOR  ON  CRT 

HALT  UNTIL  KEYBOARD   DEPRESSED 

IC=  KEYBOARD  CHARACTER 

IX=  X  CURSOR  POSITION 

IY=  Y  CURSOR  POSITON 

DSK2(14,NBUF,NWD,  JR) 

14/15  =  READ/WRITE   MASTER  FILE 

NBUF  =  BUFFER  NAME 

NWD=   NUMBER  OF  WORDS 

JR=    RELATIVE  SECTOR  OF  MASTER  FILE 

ERASE  CLEARS  CRT  DISPLAY 

EXEC(lOtNM)    CAlL  PROGRAM  STORED  IN  DISC  FILE  NM 

EX  EC ( 14, IDCON,NBUF,NWD,NAME,NDX) 

EXECQ5,  IDCON,!\IBUF,NWD,NAME,NOX) 

14/15=  READ/WRITE  DISC  FILE 
IDCON=  DISC  CONTROL  WORD 
NBUF=  START  LOCATION  OF  BUFFER 
NWD=  NUMBER  OF  WORDS  IN  BUFFER 
NAME=  NAME  OF  FILE  TO  READ/WRITE 
NDX=  STARTING  SECTOR  OF  FILE 

JFIND(  1ST,  I  END,  IVAL,  I  FID 

CONDUCTS  BINARY  SEARCH  FROM 
IFIL( 1ST)  TO  IFIL( IEND) 
RETURNS  INDEX  OF  LAST  VALUE 
LESS  THAN  OR  EQUAL  TO  IVAL 

JPLOT(K, IX, IY) 

PLOT    ON    CRT 

K    =    CONTROL    OF    POINT    OR    LINE    PLOT 

IX=  X  POSITION 

IY=  Y  POSITION 

POSAC(K)       POSITIONS  ALPHANUMERIC  CURSOR 
K=  LINE  NUMBER 

SL2RS( JMFO,NS, JS, JR) 

FIND    RELATIVE    SECTOR    FROM    MASTER 

FILE    DIRECTORY    LOCATION 

JMFO=    MASTER    FILE     INFORMATION    BUFFER 

NS=  NUMBER    OF    SIGNALS    TO    FIND 

JS=  DIRECTORY    LOCATION 

JR=  RELATIVE    SECTOR    IN    MASTER    FILE 
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INTERFACE  PROGRAM 


FTN 

PROGRAM  FACE 

DIMENSION  IFIL(3) 

DIMENSION  IFL(3) 

DIMENSION  NM(3) ,NAME(3) 

DIMENSION  NAM(3) 

COMMON  M,MOX, ISIG( 126) 

COMMON  RATE,N, IXI256) 

CCMMON  IC(30) 

COMMON  IDR(768) 

EQUIVALENCE  ( T I  MB , IC ( 7 ) )  ,  ( TT , I C (  12) ) 

EQUIVALENCE  (TOT, IC( 14) ) , ( RST, IC ( 16) ) 

EQUIVALENCE  (DIS,IC(20)) 
C 

C      THIS  PROGRAM   REFORMATS  THE  DISC  FILE  PARF  TO 
C      A  FORM  ACCEPTABLE  BY  THE   ESL  PROGRAM  RASF 
C      IT  SETS  UP  A  FILE   PMAP  TO  PASS  INFORMATION  TO  THE 
C      RASTER  DISPLAY  PROGRAM 
C 

C      VARIABLE  INFORMATION 

C         IFL    NAME  OF  FILE  CONTAINING  FORMATTED  INFO 
C         NAM    NAME  OF  FILE  TO  BE  REFORMATTED 
C         NM     NAME  OF  NEXT  PROGRAM 
C         NAME   NAME  OF  CURRENT  PROGRAM 

C         RATE, N, IX    BUFFER  FOR  SIGMAL  ENTRY  FROM   PARF 
C  RASTER  PERIOD, NUMBER  OF  ENTRIES, 

C  DELTA  TRANSITION  TIMES 

C         IC(1)    NUMBER  OF  WORDS  OF  FILE  INFO 
C         IC(2-4)    NAME  OF  FILE  WHERE  TIME  DATA  STORED 
C         IC(5)     #  SECTORS  IN  FILE  PMAP 
C  IC<6)    #  SIGNALS  IN  FILE  PMAP 

C         IC(7-8)    QUANTIZATION  INTERVAL 
C         IC ( 9-1 1 )    FILE  TYPE  CODING 
C         IC( 12-13)    SIGNAL  DURATION 

C         10(14-15)    NUMBER  OF  TRANSITION  TIMES  IN  FILE 
C         IC(i6-17)    NOMINAL  RASTER  PERIOD 
C         IDR    DIRECTORY  BUFFER  FOR  PARF 
C 

DATA  IFL(  1), IFL( 2) , IFL(3) 
X     /50115B, 405208, 2004GB/ 
C  PM     AP 

C 

DATA  NAM( 1) ,NAM( 2) ,NAM(3) 
X    /501018,51106B, 20040B/ 
C  PA      RF 

C 

DATA  NM(1 ) ,NM(2) ,NM(3 ) 
X     /51 101B,51506B,20040B/ 
C  RA      SF 


C 

C 
C 


DATA  NAME! 1) ,NAME( 2) ,NAME(3) 
X    /43101B,   41505B, 200403/ 
C  FA         CE 

CALL  EXEC  ( 14,102B,IDR,766,NAM,0) 

1  CONTINUE 

READ(5,200)   I SG 

200  FORMAT( 14) 
IF(ISG)  1,999,201 

201  CONTINUE 
DO  10  1=1.768 
IF  (  IDR( ii-ISG)  10,5,10 
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5  NDX=I*2+4 

GO  TO  20 
10  CONTINUE 
C 

WRITE  (6,101)  ISG 
101  FORMAT  ("SIGNAL  NO.  "fI5»"N0T  PRESENT") 

GO  TO  1 
C 
C 

20  CONTINUE 

CALL  EXEC  ( 14,102B,RATE,256,NAM,NDX) 

CALL  EXEC  (  15, 1028, IX  ,256,IFL,0) 
C 

C      CALCULATE  SIGNAL  DURATION 
C 

TT  =  0. 

DO  45  1=1  ,N 

IF    ( IX ( I ) )    45,45,44 
•      44    CONTINUE 

TT=TT+FLOAT( IX  (  I  )  ) 
45    CONTINUE 
C 

C  FCRMAT    FILE    INFORMATION 

C 

IC( 1)=21 

DO    40    1=1,3 

K=I+1 

IC(K)=IFL( I) 
40    CONTINUE 

IC(5)=2 

IC(6)=1 

IC(9)=52117B 

IC(10)=0 

IC( 11)=10028 

IC(18)=0 

IC(19)=ISG 

RST*-RATE 

TIMB=.0001 

RST=RST*TIM3 

TCT=FLOAT(N) 

DIS=-1. 
C 
C      SET  UP  FILE  PMAP  USING   ESL  ROUTINE 


C 


CALL  USENT  ( IC,NAME) 
CALL  EXEC  (10, NM) 
999  STOP 
END 
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MEASUREMENT  PROGRAM 


FTN 

PROGRAM  FOURD 

DIMENSION  FTX(256) 

DIMENSION  NAMRC3) 

DIMENSION    NPARO) 

COMMON  10(2048) 

COMMON  BP(256) 

COMMON  SIGU28) 

COMMON  IDMOJ  , SNUM, IDMY ( 1 1 ) 

COMMON  NAME(3) 

CCMMON  TT,T1,T2 

COMMON  RATE, PER, TB, MODE 

COMMON  IDMM(9) 

COMMON  JX,JY 

CCMMCN  MODI 

COMMON  ISG 

EQUIVALENCE  ( FTX (  I  )  ,  I D ( 1025 ) ) 

EQUIVALENCE  ( NDX , BP ( 256 ) ) 
C 

C      THIS  PROGRAM  IS  A  DRIVING  ROUTINE  WHICH 
C      CALLS  THE  MEASUREMENT  AND  INTERPOLATION 
C      SUBROUTINE  ,FAST  FOURIER  TRANSFORM, AND 
C      THE  FOURIER  SPECTRUM  DISPLAY 
C 

C      THE  FOLLOWING  VARIABLES  ARE  USED: 

C         FTX    COMPLEX  INPUT  TO  FAST  FOURIER  TRANSFORM 
C         SIG    MEASURED  SIGNAL  PARAMETERS 
C 

DATA  NAMR( 1) ,NAMR(2) ,NAMR(3) 
X      M3101B,41505B,20040B/ 
C  FA       CE 

DATA  NPAR( 1) ,NPAR(2) ,NPAR(3) 
X      /50101B,51040B,20040B/ 
C  PA      R 

C 

C      ZERO  PARAMETER   MEASUREMENTS 
C 

DO  5  1=1, 128 
5  SIG(I)  =  0. 

CALL  RESDU 
C 

C      PAUSE  FOR  RASTER  DISPLAY 
C      PRESS  Q  TO  QUIT 


C 


C 


CALL    CURSR( IC, IX, IY) 
IFUC-121B)     10,20,10 
10    CONTINUE 


C 

C      CALCULATE  128  COMPLEX  POINT  FAST  FOURIER  TRANSFORM 
C 

CALL  F0UR2(FTX, 128,1 ,-1 ) 

CONTINUE 

CALL  FORD 
C 

C      PAUSE  FOR  FOURIER  DISPLAY 
C 

CALL  CURSR( IC, IX, IY) 

IFUC-121B)  30,20,30 
C 
C      WRITE  PARAMETERS  TO  DISC  FILE  PAR 


30  CALL  EXEC(  15, 102B  ,S I G,256 , NPAR , NDX) 
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C      CALL  PROGRAM  FACE  FROM  DISC 
C 

20  CONTINUE 

CALL  EXEC  (IOiNAMR) 

ENDS 


SUBROUTINE  RESDU 

DIMENSION  IT(256) 

DIMENSION  TX(2) 

DIMENSION  FTX(256)  ,XF<2,4) ,B(5) 

DIMENSION    X( 2,123} 

COMMON  ID(2048) 

CCMMON  BP(256) 

COMMON  SIG(128I 

COMMON  IDM(3),SNUM,IDMY(11) 

COMMON  NAME(3) 

COMMON  TT,T1,T2 

CCMMON  RATE, PER, TB, MODE 

COMMON  LCT,SCX,SCR,STX,STR 

CCMMON  JX,JY 

CCMMON  MODI 

COMMON  ISG 

EQUIVALENCE  <X( 1 , 1 ), ID( 1537) ),(FTX(1),ID( 1025)) 

ECUIVALENCE( IT,ID{513) ) 

EQUIVALENCES, Tl) 

POLY(F) =(SIG(5)+F*(SIG(6)+F*(SIG(7) 
l+F*(SIG(8)+F*SIG<9)) ) ) )/PER 
C 

C  THIS  SUBROUTINE  INTERPOLATES  BETWEEN  PTS  OF  THE 

C  SCAN,  PROVIDES  MEAN  SQUARE  FITS  FOR  EACH  BIAS 

C  AND  OBTAINS  A  RASTER  VARIANCE  MEASURE 

C 

C      THE  FOLLOWING  VARIABLES  APPEAR  IN  COMMON 
C         ID(l-256)    DELTA  TRANSITION  TIMES 
C         IT       TRANSITION  TIMES 
C         X(1,K)    RASTER  PHASE  POINTS 

X(2,K)    FLOATING  PT  TRANSITION  TIME 
SIG    SIGNAL  PARAMETERS  MEASURED 
C  SNUM    TOTAL  NUMBER  OF  TRANSITIONS 

C         TT    SIGNAL  DURATION 

C         RATE f PER f TB, MODE, LCTf SCX, SCR tSTX,STR,JXfJY 
C  CONTROL  VARIABLES  FOR  DISPLAY 

C   .       ISG    SIGNAL  IDENTIFIER 

C         XF(2,4)   4  POINTS  USED  FOR  CUBIC  INTERPOLATION 
C         B    COEFFICIENTS  OF  INTERPOLATED  FIT 
C 

C      INITIALIZATION 
C 

NTOP=IFIX(SNUM) 

LT  =  0 

RX=JX 

RY=JY 

LX=900 

LY=500 

SCAX=FLOAT(LX)/TT 

RAT2=PER/2. 

RAT4=RAT2/2. 
C 

C  CONVERT    DELTA    TIME    TO    TIME 

C 

DO    200    I=1,NT0P 

IF(ID( I ))  200,200,201 
201  IDEL=ID(I) 

LT=LT+IDEL 
200  IT(I)=LT 
C 

C      FIND  STARTING  INDEX  FOR  MINIMAL  PHASE  DIFFERENCE 
C 

JST  =  JFIND( 1, NT  OP, 400, I T) 
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SIG(2)=FL0AT( I T ( NTOP-2) -400 ) 
DELT=  SIG(2)/128. 
SIG(2)=1./(SIG(2)*TB) 

N0X=1 

NSWCH=1 

OHSQ=0. 

JTOP=0 

KST=1 
999  CONTINUE 

PHM=0. 

PHSQ=0. 
C 

C      INTIALIZE  RASTER  PHASE  POINT 
C 

J  =  JST 
210  J  =  J-2 

IF(J-l)  212,212,209 

212  PHH=0. 
TST=0. 

GO  TO  214 
209  IF(IDU))  210,210,213 

213  TST=IT< J) 
PHH=AMOD(TST,PER) 
IF(PHH-RAT2)  214,215 

215  PHH=PHH-PER 

214  X(lfl)=TST 
X(2,1J=PHH 
DELL=0. 
PHM=PHM+PHH 
PHSQ=PHSQ>PHH*PHH 
N=TST/PER 

J  =  l 

PhL=PHH 
C 

C      CALCULATE  MINIMAL  PHASE  DIFFERENCE  POINTS 
C 

DO  300  I=JST,NTOP, 2 

IF(IDU))  300,300,302 
302  TST=IT( I) 

41  N=N+1 
PHH=TST-FLOAT(N)*PER 
DEL=PHH-PHL 
IF(ABS(DEL)-RAT2)  42,42,41 

42  N  =  N-1 
C 

C  UPDATE    MEANS    AND    VARIANCE    ESTIMATES 

C 

PHK=PHM+PHH 

PHSC=PHSQ+PHH*PHH 

OHSQ=OHSQ+DEl*DEL 

IF(ABS(DEL)-RAT4i     305,707 
707    DEL=DELL*.5 

PHH=PHL+DEL 
305  J=J+1 

X( 1, J)=TST 

X(2, J)=PHH 

DELL=DEL 

PHL=PHH 
300    CONTINUE 
C 

c 

C      CALCULATE  MEAN  AND  VARIANCE  FOR  BIAS  CONDITION 
C 

JTOP=JTCP+J 

JP  =  J 

XNUM=JP 

XNUM1=  JP-1 

SIG(35)  =  SIG(35)+(PHSQ-PHM*PHM/XNUMJ/XNUM1 

PHM1=PHM/XNUM 
C 
C 
C      CALCULATE  POLYNOMIAL  MEAN  SQUARE  FIT 
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CALL  CURV(JP) 
C 

C      INTERPOLATE  BETWEEN  POINTS 
C 

K=l 

CALL    CURFCXdtKJ  ,3,B) 

XP=400. 

DC    601    1=1,128 
607     IFUP-X  (  l,K+2))     602»603 

602  CONTINUE 
DX=XP-X( 1,K) 

FTX(KST)=(B( 1 )+DX=MB(2)+DX*(B(3)+DX*B(4)  )) )/PER 

IX=XP*SCAX+RX 

IY=(FTX(KST)+.5)*FL0AT(LY/+RY 
C 

C  DISPLAY    INTERPOLATION 

C 

CALL    JPLOT(0,IX-4,  IY-6) 

CALL    ALPHA 

CALL    CHOUT(52B) 
C 

C  INPUT    FOR    FAST    FOURIER    TRANSFORM 

C 

FTX(KST)=FTX(KST )-POLY(XP) 

KST=KST+2 

XP=XP+DELT 

GO    TO    601 

603  IF(K+4-JP)     604,602,602 

604  CONTINUE 
K  =  K  +  1 

CALL    CURF(X( If  K)  ,3,B) 

GO    TO    607 
601    CONTINUE 

CALL    CH0UT(-1 ) 
C 

C      PLOT  POLYNOMIAL  MEAN  SQUARE  FIT  DATA 
C 

XP=400. 

DO  20  1=1,128 

IX=XP*SCAX+RX 

IY=(POLY(XP)  +.5)*FL0AT(LY)+RY 

CALL    JPLGT(0,IX-4, IY-6) 

CALL    ALPHA 

CALL    CHOUT(53BI 

XP=XP+DELT 
20    CONTINUE 

CALL    CHOUT(-l) 
C 

C   CHECK  TO  SEE  IF  BOTH  BIAS  CONDITIONS  ARE  PLOTTED 
C 

NSWCH=-NSWCH 

IF(NSWCH)  700,701 

700  CONTINUE 
XP=400. 
PHM2=PHM1 
JST=JST+1 
KST  =  2 

SIG<4)=SIG(3) 
DO    702     1=1,5 
Kl=I+4 
K2=I+9 

SIG(K2 )=SIG(K1) 
702    CONTINUE 
GO    TO    999 

701  CONTINUE 
C 

C      CALCULATE  MEAN  BIAS  AND  INTRINSIC  VARIANCE 
C 

XNM=JT0P*2 

SIG(20)=  ABS(PHM2-PHM1) 

SIG(36)=  OHSQ/XNM 
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c 
c 

C      CALCULATE  LAGUEKRE   COEFFICIENTS 


C 


P  =  80. 

TEND=800. 

CALL  LAGER( I D , PE R, P, TEND , S I G( 15) ) 

RETURN 

END 


SUBROUTINE  CURV(NPTS) 

DIMENSION  X(2,128) 

DIMENSION  B(5) 

DIMENSION  A(5,6) 

COMMON  ID12048) 

COMMON  BP(256) 

COMMON  SIG(128) 

COMMON  IDMMYU6) 

COMMON  NAME(3) 

COMMON  TT,T1,T2 

COMMON  RATE, PER, T6, MODE 

COMMON  LCT,SCX,SCR,STX.,STR 

COMMON  JX  JY 

CCMMON  MODl,ISG 

EQUIVALENCE  ( X( 1 , 1 ) , ID ( 1537 ) ) , ( B ( i ) , SI G ( 5) ) 
C 

C      THIS  SUBROUTINE  CALCULATES  THE  COEFFICIENTS 
C      OF  A  POLYNOMIAL  MEAN  SQUARE  FIT  OF  ORDER  LESS 
C      THAN  FOUR, USING  UP  TO  128  SAMPLE  POINTS 
C 

C  X(2f128J    SAMPLE  POINTS 

C  SIGC5-10)    POLYNOMIAL  COEFFICIENTS 

C 
C      SPECIFY  ORDER  OF  POLYNOMIAL  TO  BE  FIT 


C 


NUM=2 

DO  20  J  =  l,6 
DO  20  1=1,5 
A(I,  J)=  0. 
20  CONT INUE 

A(1,2)=FL0AT(NPTS) 

DO  100  I=1,NPTS 

YT=  X(2,I) 

X1=X( 1,1) 

X2=  X1*X1 

X3=  X2*X1 

X4=  X3*X1 

X5=  X4*X1 

X6=  X5*X1 

X7=  X6*X1 

X8  =  X7*X1 

A<1,1)=  A(  1, 1)+YT 
A<  1,3)=  A( 1,3)+X1 
A(l,4) =A( 1,4)*X2 
A(1,5)=A( 1,5)+X3 
A(1,6)=A( 1,6)+X4 
A(2,1)=A( 2,1 )+Xl*YT 
A(2,6)=A( 2,6)+X5 
A(3,1)  =  A(3,1  )+X2*YT 
A(3,6)  =  A(3,6)«-X6 
A(4,1)=A(4,1  )*X3*YT 
A(4,6)=A(4,6 )+X7 
A(5,1)=A( 5,1)+X4*YT 
A(5,6) =A(5,o)+X8 
100  CONTINUE 

DO  25  1=2,5 
DO  25  J=2,5 
A(I,J)=A( I-1,J*1) 
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c 
c 


,LMT 
,LCT 
I,K)/A( 


I  ,M) 


7,45,45 


25  CONT  INUE 
LCT=  NUM+2 
LMT  =  NUM+1 
DO  35  1=1, NUM 
DO  35  J  =  l,  I 
L  =  J+l 
K=I+1-J 
M  =  K  +  1 

FAC=  A(K,L)/A(M,L) 
CALL  ROW(A,NUM,M,FAC»K) 
CONTINUE 
M=  LCT 
DO  45  1=1 
DO  36  K=l 
A(I ,K)=A< 
CONTINUE 
IF( I-LhT) 
CONTINUE 
DO  38  J=I 
MK=J+1 
FAC=A( MK, 
CALL  ROW( 
CONTINUE 
M  =  M-l 
CONTINUE 
J  =  LCT 

DO  46  1=1, 
J  =  J-1 

B(J)  =  A(  I, 
CONTINUE 
IFCLCT-5) 
DO  47  I=LC 
B( I)=  0. 
CONTINUE 
RETURN 
END 


35 


36 
37 


38 
45 


46 

70 

47 
75 


,NUM 

M) 

A,  NUM 


LMT 


I ,FAC,MK) 


1) 

70,70,75 

T,5 


SUBROUTINE  ROW  ( A, NUM , IRM, FAC» IRS ) 
DIMENSION  A (5, 6) 
DIMENSION  TEMP (b) 
KJ=NUM+2 
DO  50  1=1, KJ 
TEMP(I)=  FAC*A(IRM,I) 
50  CONTINUE 

DO    100     1=1, KJ 

A(  IRS, I  )=    A(  IRS,  D-TEMPCI  ) 
100    CONTINUE 
RETURN 
END 
ENDS 


C 

c 
c 
c 
c 


c 
c 
c 
c 


SUBROUTINE  L AGE R ( I D, RATE , P, TEND  SL) 
DIMENSION  10(256), SL(5) 

THIS  SUBROUTINE  CALCULATES  THE  LAGUERRE  POLYNOMIAL 
COEFFICIENTS  FOR  THE  FIRST  FIVE  POLYNOMIALS  USING 
TRAPEZOIDAL  INTEGRATION 

DO  5  1=  1,5 
SL(1)=0. 

CALCULATE  MINIMAL  PHASE  DIFFERENCE  RASTER  POINTS 
USING  BOTH  BIAS  CONDITIONS 

PHL=0. 
TT=  0. 
N=  0 
P=  l./P 
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RAT2  =RATE/2. 
DO  10  1=1,256 
IF(  ID(  I  ) )  10,  10,  11 

11  TT=TT  +  FLOAT( ID( I  )  ) 
IF(TT-TEND)  12,12,13 

12  DEL  =  FLOAT(  ID( I )  ) 

14  N=N+1 
PHH=(TT-FLOAT(Ni#RATE) 
OIF=  PHH-PHL 

IF(ABS( DIF)-RAT2)  15,15,14 

15  N  =  N-1 
C 

C  CALCULATE    CONTRIBUTIONS    TO    LAGUERRE    POLYNOMIALS 

C 

K=l 

T1=TT-DEL 

30  K=K-1 
DIF  =  PHL 
RT-P*T1 

REXP=.5*DEL*EXP( -RT) 
SL(1)=    SL(  1)+REXP*DIF 
SL(2)=SL(2)+(2.*RT-1.)*REXP*DIF 
SL(3)=SL<3)+( l.+RT*(-4.+RT*(2.) ) )*REXP*DIF 
SL(4)=SL(4)  f  <-1.4-RT*<6.+RT*(-6.  +RT*( 4. / 3 . ) ) )  )*R£XP*DIF 
SL(5)=SL(5)+<1 .+RT*{-8.+RT*(12.+RT*(-16./3.+RT*2./3. )) 

1  ))*REXP*DIF 

1    *REXP*DIF 

T1  =  TT 

PHL=PHH 

IF(K)     31,30,30 

31  CONTINUE 
10    CONTINUE 

C 

C  NORMALIZE 

C 

13  PN0RM=SQRT(2./P) 

DO    20     1=1,5 
20    SL(I)=SL( I)*PNORM 
RETURN 
END 


SUBROUTINE    FORD 

DIMENSION    FX(256) 

DIMENSION    NM(3) 

DIMENSION    FTX(2  56) 

COMMON  ID (2048) 

COMMON  BP(256) 

COMMON  SIG(128) 

COMMON  IDM(3) ,SNUM, IDMY(ll) 

COMMSN  NAME(3) 

COMMON  TT,T1,T2 

COMMON  RATE, PER, TB, MODE 

COMMON  IDMM19) 

COMMON  JX,JY 

CCMMCN  MODI 

COMMON  ISG 

EQUIVALENCE  ( FX (  1 ) ,  ID( 1 537 )  ) 

EQUIVALENCE  ( FTX ( 1 ) , I D ( 1025 ) ) 

EQUIVALENCE  ( NDX , B P ( 256 ) } 
C 

C      THIS  SUBROUTINE  CALCULATES  AND  PLOTS  FOURIER 
C      MAGNITUDES  FOR  THE  64  TERMS 


C 


RX  =  JX 
RY  =  JY 
LX=  900 
LY  =  500 
XMAX=  200. 
YMAX=4. 
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XMX4=  XMAX/4. 
YMX2  =  YMAX/2, 
SIG(1)=  RATE 
STX=0. 
STY=0. 

SCX  =  FLGAT(LX)/XMAX 
SCY=  FLOAT(LY)/YMAX 
CALL  ERASE 
CALL  P0SAC(2) 

WRITE(6,100)  ISG,SIG<37) ,SIG(2) 
100  FGRMATt "FREW  MAG  PLOT  FOR  SIG  NO. ",15/ 
1, "VARIANCE  OF  RE S I DU£  =  " , F 10  .6 . 
2  "   RESOLUTION=  ",F10.6) 

CALL  AXIS< JX, JY,SCX,SCY,STX,STY, 
1  XMX4,YMX2,XMAX,YMAX) 

C 
C      CALCULATE  FIRST  64   MAGNITUDE  COMPONENTS 


C 


DO  40  1=1,63 
J=2*I+1 
K=257-2*I 
L  =  I  +  1 
JI=J+1 
KI=K+1 

FX(L)  =  (  SQRT(FTX( J ) **2+FTX( K ) **2  +  FTX( J  I ) **2 
1  +FTX(KI  )**2) )/2. 
40  CONTINUE 

FX(1)=( SQRT(2.*(FTX(1)**2+FTX(  2)**2) ) )/2. 


c 
c 

c 

c 

PLOT    SPECTRUM 

50 

DO    50    1=1,64 

IX=FLOAT( I~1)*SIG(2J*SCX+RX 

I Y=FX(  I  )*SCY+RY 

CALL    JPLGTU-I  ,  I  X  ,  I  YJ 

CONTINUE 
CALL    ALPHA 
CALL    CH0UT(-1 ) 

c 

10 

DO    10    1=1,64 
K=I+39 
SIG(K)=FX( I) 

CONT INUE 

RETURN 

END 
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STATISTICS  AND  SCALING  PROGRAM 


FTN 

PROGRAM  CSTAT 

DIMENSION  SPAR(50) , JW( 50) 

01  MENS  ION    NPAR  (3)  ,NCLS(3  ),  NL  IB  (3) 

COMMON  JSAVE(50) 

COMMON  JMFO  (128)  ,JMFD(768),  MFE(334) 

COMMON  IPD(256) ,SIG( 128) 

COMMON     ICDU28) , ISG( 128) 

COMMON    JBUF(4) iCLMN(50)fCLVA(50), JNPAC50),     JPNUM(50), 
1  JCDO(80) 

COMMON  JGL06<4) , GM( 50) , GV( 5 0) , G TV( 50 ) , J8ND( 50 ) , 
1        JGDM130) 

EQUIVALENCE  ( SPA R( 1 ) , MFE ( 5 )  ) , ( JW (  1 ) , MFE (  105 ) ) 

EQUIVALENCE  ( IDMU, JBUF( 1)  )  , ( NEN T, JBUF ( 3 ) ) 
C 

C      THIS  PROGRAM  PLACES  MEASURED  SIGNAL  PARAMETERS  IN  THE 
C      SIGNAL  FILES  OF  THE  MASTER  FILE.   CLASS  STATISTICS 
C      ARE  CALCULATED  AND  STORED  IN  THE   MASTER  FILE  CLASS 
C      ENTRIES. 
C 
C      THE  FOLLOWING  VARIABLES  ARE  USED: 

JSAVE    STORAGE  FOR  OLD  PARAMETERS 
C         JMFO    INFORMATION  ON  MASTER  FILE  STRUCTURE 
C  JMFD    MASTER  FILE  DIRECTORY 

C         MFE    BUFFER  FOR  MASTER  FILE  ENTRY 

IPD    PAR  FILE  DIRECTORY 
C         SIG    BUFFER  FOR  PAR  FILE  ENTRY 
C         ICD    CSIGF  FILE  DIRECTORY 

C         ISG    BUFFER  FOR   SET  OF   SIGNAL  NUMBERS  OF  SAME 
C  CLASS 

C         JBUF-JCDO    BUFFER  FOR  ENTRY  TO  CLASS  FILE 
C         JGLOB-JGDM    BUFFER  FOR   ENTRY  TO   GLOBAL  FILE 
C 
C 

DATA  NPAR(l) ,NP AR < 2 ) , NPAR ( 3 ) 
1    /50101B,  51040B,  200408/ 
C  PAR 

DATA  NCLS( 1) ,NClS(2) ,NCLS(3) 
1    /41523B,  44507B»   43040B/ 
C  CS        IG      F 

DATA  NLIB( 1) ,NL  I  B { 2 ) , NL I B  (  3 ) 
1   /46501B, 51524b, 430408/ 
C  MA       ST      F 

C 

C       READ  IN  DIRECTORIES  FOR  FILE   PAR,CS I GF , MASTF 
C 

CALL  EXEC( 14, 1028, IPD  ,2  56 , NPAR , 0 ) 

CALL  EXECt  14, 1028,  ICD  ,  1 26 , NCLS , 0 ) 

CALL  EXEC( 14,10  28, JMFO, 12  3, NL IB, 0) 

NWD=JMF0(2) 

CALL  DSK2(  14,JMFD,NWD, 1 ) 

CALL  SL2RS(JMF0,1 , 1,LG) 

CALL  DSK2( 14, JGLOB , 384 , LG ) 
C 
C         ZERO  GLOBAL  MEANS  AND  STND  DEV 


C 


DO  10  1=1,50 
GM(I )=0. 
GV( I )=0. 
GTV( I)=0. 
10  CONTINUE 
JGL08( 3)=0 
CALL  ERASE 
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c 
c 


CALL  POSAC(O) 
WRITE (6,113) 
113  FORMAT ( "ENTER  INPUT  DEVICE  (5  FOR  CR,  6  FOR  TEK)  " ) 


READdt*]  LU 

READILU,*)  N 
C 

C   CHANGE  PARAMETER  NAMES  IN  MASTER  FILE  INFO  SECTOR 
C 

DO  90  I  1  =  1, N 

Kl  =  11+14 

JSAVE(  U)=JMF0(K1  ) 

READ(LU,*)  JMFO(Kl) 
90  CONTINUE 

JMFO( 14)=N 
1  CONTINUE 

CALL  ERASE 

CALL  POSAC(O) 

WRITE(6,100) 

100  FORMAT ("ENTER  CLASS  NO."/) 
READ(6,*)  NUM 

C 

C         READ  IN  SIGNAL   IN  CLASS   SPECIFIED 

C 

CALL  EXEC(14,102B,ISG,128,NCLS,NUM) 

IDNO=ICD(NUM) 
C 

C  FIND  CLASS  ENTRY  IN  MASTF  AND  READ  INTO  JBJF 
C 

DO  70  JM=1,NWD 

IF    (  JMFD(  JMKIDNO)     70,71,70 

70  CONTINUE 
GO    TO    1 

71  CONTINUE 

CALL    SL2RS(JMF0, 1, JM, JR) 

LOCLI=JR 

CALL  DSK2( 14, JBUF,384, LOCLI ) 
C 

C         ZERO  FILE  MEANS  AND  VARIANCE 
C 

DO  11  1=1,50 

JNPA( I) =0 

CLMN( I ) =0. 

CLVA( I)=0. 
11    CONTINUE 

NENT=0 
C 

C         COMPUTE  MEAN  +  VARIANCE  OF  PARAMETERS. 
C 

DO  20  K=l ,128 

IFUSG(K)  )  21,21,22 
C 

C      FIND   RELATIVE   SECTOR   OF   ISG 
C 

22  DO  30  J=l,256 

IF(  IPD( J)-ISG(K)  )  30,23,30 
30  CONTINUE 

WRITE(6,101)  ISG(K) 

101  FURMATt "SIGNAL  NO.", 16,"  NOT  FOUND"/) 
GO  TO  20 

23  NENT=N£NT+1 
JGLOB(3)=JGLOB(3)+l 
NDX=2*J+2 

CALL  EXEC(14,  102B , SI G, 256 , NPAR, NDX) 

CALL  TRNSGUSG(K)  ,  IDNO) 
C 

C   CALCULATE    SUM  +   SUM  OF   SQUARES 
C 

DO  40  1=1,50 

A=SPAR( I) 

GM( I )=GM( I )+A 
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GV(I )  =  GV( I )f A*A 
IF( JW( I ))  40,40,50 

50  CONTINUE 

JNPA(  I  J  =JNPA( I)  +  i 

XNUM=FLOAT( JNPA( I) ) 

XNUM1=XNUM-1 . 

XNUM2=XNUM-2. 

CLMN( I) =  (CLMN( I ) *XNUM1 +SP AR( I) J/XNUM 

IF(JNPA( I )-l)  51,51,52 

52  Continue 

clva( i ) =xnum2/xnum1*clva( i ) + ( ( sp ar ( i ) -clmn (  i )  ) **2  )  / 

1  XNUM 

51  CONTINUE 
40  CONTINUE 

20  CONTINUE 

21  CONTINUE 

CALL  CURSR( IC, IX, IY) 
DO  80  1=1,50 
JM=I+14 

JPNUMC I )  =  JMF0(  JM) 
80  CONTINUE 

CALL  0SK2( 15,J3UF,384,L0CLI> 
CALL  ERASE 
CALL  POSAC(O) 
WRITEC 6,102) IDNO 

102  FORMAT ("CLASS  NAME    "  ,16/ 

1"  PAR     MEAN      VARIANCE    fcUSED"/) 
WRITE ( 6,103)  ( I,CLMN( I ),CLVA(I) , JNPA(I) ,1=1 ,50) 

103  FORMAT* 2( I  6, 2E1 2 .4 , 14, 2X)  ) 
CALL  CURSRt IC,IX,IY) 
IFUC-121B)  1,999,1 

999  CONTINUE 
C 

C         CALCULATE  GLOBAL  MEAN  AND  STD  DEV 
C 

CALL  ERASE 

CALL  POSAC(O) 

WRITE(6,800) 

800  FORMAT ("GLOBAL  VALUES"/, 

1"PAR     GLOBAL  MEAN     STD  DEV     TY P  DEV     JBND") 

801  F0RMAT(6X,I4,3(2X,E12.4) ,1X,I4) 
C=  JGLUb(3) 

DO  802  1=1,50 

A=GM( I ) 

B=GV( I ) 

A=A/C 

D=SQRT(  <B-C*A*A)/<C-1.)) 

IF(D)  803,803,804 

803  D=l. 

804  GM( I )=A 
GV(I )=D 
GTV( I )=.2*D 
JBND( I)=2 

802  CONTINUE 
C 

C        SET  UP  VALUES  IN  SECTOR  0  OF  MASTER  FILE 
C 

JMF0(5)=-1 

JMF0(6i=5 

CALL  EXEC(15,102B, JMFO, 128 , NLI B , 0 ) 
C 
C        PRINT  OUT  VALUES  OF  GLOBAL  FILE 


C 


WRITE (6, 801) (JMFO(  1  +  14) ,GM( I ) ,GV(  I  ),GTV(  I)  , 
1JBND( I ) ,1=1,50) 
CALL  DSK2( 15 , JGLOB , 384, LG) 
STOP 
END 


149 


SLBROUTINE  TRNSG ( I S , LDUM) 

DIMENSION  WSC115) 

DIMENSIGN  SPAR(50) T JW( 50) 

DIMENSION  NPAR(3),NCLS(3),NLIB(3) 

COMMON  JSAVE(50) 

COMMON  JMFO  (123)  ,JMFD(768),  MFE<384) 

COMMON  IPDI256) ,SIG( 128) 

COMMON  ICD(L28)  T  ISGU28) 

COMMON  JBUFU)  ,CLMN(50)  ,CLVA(50)  ,  JNPA(50)  ,  JPNUM(50), 
1       JCDO(80) 

COMMON  JGL0B<4) , GM(50J , GV(50) , GTV(50) ,JBND( 50) , 
1        JGDM(30) 

EQUIVALENCE  ( SP AR ( 1 ) , MFE ( 5 )  ) , ( J W (  1 )  ,MFE (  1 05 ) ) 

EQUIVAL  ENCE  ( IDNO, JBU F ( 1 ) )  , ( NEN T , JBUFi 3 )  ) 
C 

C      THIS  SUBROUTINE  CALCULATES  PARAMETERS  AND  STORES  THEM 
C      IN  THE  MASTER  FILE  BUFFER   MFE 
C 

c 

DATA  WSC( 1) »WSC( 2) ,WSC(3) ,WSC(4) 
1    /1000. ,1. tlOOO. ,1000./ 

DATA  WSC(5) ,WSC(6) ,WSC(7) ,WSC(8) ,WSC(9) 
1    /I. , 100.,1.E06,1.E09,1.E13/ 

DATA  WSC( 10) ,WSC( 11) , WSC(12) ,  WSCU3) ,WSC(14) 
1    /l. , 100., I.c06, l.E09,l.E13/ 
C 
C 

C      THIS  SUBROUTINE  SCALES  PARAMETERS  AND  STORES  THEM 
C      IN  THE  MASTER  FILE  BUFFER   MFE 
C 
C 

NTOP=JMF0( 2) 

DO  5  1=1, NTOP 

IF  (IS-JNFD(D)  5,6,5 

5  CONTINUE 

DO  7  1=1,50 
JW(I)=0 
7  CONTINUE 

WRITE (6, 100) 
100  FORMAT( 15, IX, "NOT  IN  MASTF") 
RETURN 

6  CONTINUE 
K  =  0 

CALL  SL2RS( JMFOt 1 , I, JR) 

CALL  DSK2( 14, MFE ,384, JR) 

MFE(1)=IS 

MFE(3)=-JBUF(1) 

MFE(4)=0 

K  =  K  +  1 

JWT  =  1 
C 

C      RATE 
C 

SPAR(l) =  SIG( 1)*WSC(1) 
C 

C       RESOLUTION 
C 

SPAR(2)=  SIG(2) 
C 

C  POLYNOMIAL  MEAN  SQUARE  FIT  COEFFICIENTS 
C 

SPAR(3)=  SIG(5)*WSC(5) 

SPAR(4)=SIG(6)*WSC(6) 

SPAR(5)=  SIG( 7)*WSC( 7) 

SPAR(6)  =  SIG(  10)*WSC(  10) 

SPAR{ 7) =SIG( 11 J*WSC(11 ) 

SPAR(8)=  SiG(12)*WSC(12) 
C 

C       LAGUEkRE  POLYNOMIAL  COEFFICIENTS 
C 
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SPAR(9)=  SIG(15) 

SPAR(10)=  SIG(16) 

SPARdl  J=  SIG(  17) 

SPAR( 12 )=SIG(18) 

SPAR(13)=  SIG(19) 
C 

C       MEAN  BIAS 
C 

SPAR(14)=  SIG(2Q) 
C 

C      VARIANCE 
C 

SPAR(15)=  SIG(35) 
C 

C      INTRINSIC  VARIANCE 
C 

SPAR(16)=  SIG(36) 
C 

C  NONPARAMETRIC    CORRELATION 

C 

SPARU7)  =    l.-SIG(36)/SIG(35) 
C 

C       FOURIER  POWER  SPECTRUM  USING  LINEAR  INTERPOLATION 
C 

K=  18 

J=40 

L  =  l 

RES=SIG(2) 

SPAR(K)=ABS(SIG( J)  ) 

DO  60  1=1,32 

K=K  +  1 
63  CONTINUE 

DIFF=FLOAT( I )*2.-FL0AT( L)*RES 

IF(DIFF)  61,61,o2 

61  SPAR(K)=SIG( J+1)+(SIG(J+1)-SIG( J))*DIFF/RES 
GO  TO  60 

62  L=L+1 
J  =  J  +  1 

GO  TO  63 

60  CONTINUE 

256  CONTINUE 

DO  70  I=1T50 
JW( I)=l 
70  CONTINUE 

CALL  DSK2( 15,MFE,384, JR) 

RETURN 

END 
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